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Summary 


A detailed study has been made of face and annular seals under conditions where 
boiling, i.e., phase change of the leaking fluid, occurs within the seal. Many seals operate 
in this mode because of flashing due to pressure drop and/or heat input from frictional 
heating. We mention high pressure, water pumps, industrial chemical pumps, and 
cryogenic pumps as a few of many applications. The initial motivation for this work was 
the LOX-GOX seals for the space shuttle main engine, but the study has been expanded to 
include any face or annular seal where boiling occurs. 

We have discussed some of the distinctive behavior characteristics of two-phase 
seals, particularly their axial stability. While two-phase seals probably exhibit instability to 
disturbances of other degrees of freedom such as wobble, etc., under certain conditions, 
such analyses are too complex to be treated at present. Since an all liquid seal (with parallel 
faces) has a neutral axial stiffness curve, and is stabilized axially by convergent coning, 
other degrees of freedom stability analyses are necessary. However, the axial stability 
behavior of the two-phase seal is always a consideration no matter how well the seal is 
aligned and regardless of the speed. Hence, we might think of the axial stability as the 
primary design consideration for two-phase seals and indeed the stability behavior under 
sub-cooling variations probably overshadows other concerns. The main thrust of this 
work has been the dynamic analysis of axial motion of two-phase face seals, principally the 
determination of axial stiffness, and the steady behavior of two-phase annular seals. 

The main conclusions are that seals with two-phase flow may be unstable if 
improperly balanced. Detailed theoretical analyses of low (laminar) and high (turbulent) 
leakage seals are presented along with computer codes, parametric studies, and in particular 
a simplified PC based code that allows for rapid performance prediction: calculations of 
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stiffness coefficients, temperature and pressure distributions, and leakage rates for parallel 
and coned face seals. 

A simplified combined computer code for the performance prediction over the 
laminar and turbulent ranges of a two-phase seal is described and documented. 

This report summarizes the analyses, results, and computer codes, but for more 
details the reader is referred to the more complete detailed studies presented in the various 
papers and reports listed in Chapter 8. 
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opening force [N] 
nondimensional opening force 

closing force [N] 

mass velocity, pu, [kg/m 2 s] 

film thickness (seal clearance) [m] 

nondimensional film thickness 

film coefficient for heat transfer [J/(s.m 2 K)] 

specific enthalpy [J/kg] 

thermal conductivity [W/(m.K)], inlet loss coefficient 

mass leakage rate [kg/s]; 

nondimensional mass leakage rate 

Mach Number 

pressure [Pa] 

nondimensional pressure 

torque [Nm] 

rate of heat conduction per unit area into the seal plates from the 
fluid [J/(s.m 2 )] 
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shaft radius 

ideal gas constant [J/(kg.K)] 
temperature [K] 
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specific volume [m^/kg] 
relative axial velocity of the seal rings [m/s] 
velocity in 6 - direction [m/s] 
control volume shaft work 
nondimensional radial distance 
axial coordinate [m] 
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CHAPTER 1 
INTRODUCTION 

Seals are mechanical devices used to restrict leakage of fluids, for example, when a 
rotating shaft penetrates a stationary housing which encloses the pressurized fluid. The 
tolerable leakage rate depends on the nature of the sealed fluid; leakage of expensive, toxic, 
corrosive, explosive or flammable fluids must be reduced to a minimum. The life and 
reliability of seals are also of major concerns amount the users to reduce equipment and 
process downtime. Sometimes when system redundancy is kept at a bare minimum, for 
example in airborne and space vehicles, a seal failure could cause serious system 

malfunctioning. 

In January 1986, the whole world suddenly became aware of the crucial importance 
of fluid sealing technology when the US shuttle "Challenger" tragically exploded shortly 
after leaving the launch pad. A joint sealed by rubber O-rings had failed. This episode had 
the characteristics of many a sealing problem. The component involved was of relatively 
low-value in its own merit, but the consequential cost of failure was totally 
disproportionate. The failed O-ring was a static seal and much less complicated in 
operation than the dynamic seals discussed here. 

The Figure 1-1 [1] gives an overview of different types of industrial sealing devices 
available. Among the different kinds of seals, 'Mechanical End Face Seals (also simply 
called Face Seal) are the dominating category of major industrial seals and have been given 
special and extensive considerations. Fluids that need to be sealed range from water, 
petroleum products, oil, natural gas, air and toxic chemicals to cryogenic fluids like liquid 
oxygen and hydrogen (Space Shuttle Turbo Pumps). These seals may handle pressure up 
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to 5000 psi (~350 atm.), temperatures up to 1000°C and a rotational speed of up to 60,000 
RPM. 

A good treatise on mechanical face seal designs, basic configurations, operation and 
lubrication mechanisms has been given by Ludwig & Griener [2, 3]. Figure 1-2 shows the 
schematic diagram of a face seal. The primary sealing is accomplished by a nonrotating 
ring (called primary seal ring or stator) that bears against the face of a rotating ring (called 
seal seat or rotor) mounted on the shaft Occasionally co-rotating and counter-rotating seals 
(advanced aircraft engines) are encountered where both the rings are rotating. Between the 
stator and the housing, there are multiple springs which give it the flexibility in the axial 
and two angular modes about orthogonal diametrical directions. Secondary seals are 
provided between the stator ring and the housing. Typically these seals are elastomeric O- 
rings. They self-energize under pressure and tend to fill in the asperities and voids on the 
surfaces in contact and hence minimize leakage through secondary sealing surfaces. 
Successful operation of seals requires satisfaction of seemingly competing demands. In 
order to reduce wear and maintain integrity of the sealing surfaces, it is desirable, if not 
essential, to achieve and maintain separation of faces by a lubricating film. At the same 
time face separation must be kept extremely small (~ 2-3 fi m) in order to minimize leakage. 
These requirements must be dynamically met in changing operating conditions and in the 
presence of machinery vibrations. 

Figure 1-3 shows the typical forces a stator experiences. This is an 'outside- 
pressurized' arrangement with the high pressure fluid at the Seal OD. This configuration 
offers a few advantages as the centrifugal forces tend to retard leakages and centrifuge solid 
particles away from the sealing surfaces giving a self-cleaning feature. The sealed fluid 
leaks through the gap between the seal rings and pressure drops due to friction and inertia. 
A typical pressure profile, P (r ) , for an axisymmetric gap is shown in the figure. If the 
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gap is nonaxisymmetric, hydrodynamic pressure will build up inside the seal due to relative 
tangential motion of the seal faces and the pressure profile becomes also a function of the 
circumferential location, i.e. P (r . 6 ). For an axisymmetric gap, there is no 
hydrodynamic pressure generation and the pressure distribution inside the seal is same as if 
both the seal faces are stationary (if centrifugal inertia effects are neglected). The pressure 
distribution, so obtained, often referred to as hydrostatic component. This fluid pressure 
P, tends to open the seal gap. On the other hand, the axial loading from the sealed fluid 
pressure and the spring force, F s , acts behind the stator and tends to close the seal. The 
expressions for the opening and closing forces are given below. 

F Q = f 2n f'°P{ r,0)rdrd0 

F c = Jt( r o - r bal) P o +F s + rc^bal - r ?) P i 

'r ba l ’ is called the ’balance radius’ by which the closing force, F c , can be controlled. If 
the closing force, F c , is equal to the opening force, F a , at the operating point, then a 
lubricating fluid film is maintained at the seal interface. Under this situation, the seal 
operates in a ’non-contacting mode’ and is called a 'balanced seal.' The corresponding 
clearance is called the 'operating clearance or film thickness.’ On the other hand, if the 
closing force is greater than the available opening force, the asperity contact takes place at 
the seal interface and the force balance is achieved with the help of the mechanical contact 
pressure. In this case, the seal operates in a 'contacting mode' and is called an 'unbalanced 
seal.' The contacting seals are supposed to operate with a minimal contact pressure, 
otherwise, heavy wear at the surfaces would cause premature seal failure. These seals are 
generally used for low to moderate pressure services and noncontacting seals for high 
pressure applications. In chemical and petrochemical industries, the contacting mode is 
primarily chosen for almost all sealing applications to reduce leakage of hazardous fluid as 
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much as possible even at the expense of seal life. The third situation arises when the 
opening force exceeds the applied closing force. In that case, the seal pops open causing 
high leakage and seal failure is said to have taken place. For a given design with a certain 
balance radius, rbal, the closing force is constant for a given operating pressure, whereas 
the opening force is dependent on the gap geometry and speed. The information most 
useful to the seal designers is the ’Opening Force vs. Nominal Clearance' curves, typically 
known as ’F-h’ curves, for different speeds and system pressures. One necessary 
requirement for a stable and successful seal operation is to have a negatively sloped ’F-h’ 
curve around the operating point (which means positive film stiffness); otherwise, seal 
faces will collapse and give unacceptable contact load and rapid wear. Examples of typical 
’F-h' curves are shown in Figures 1-9 and 1-10 and they are explained later. 

The seal lubricating film is usually very thin (in the range of few microns) and, 
therefore, very small irregularities, thermal and pressure distortions, and face runout 
motions can have a dramatic effect on seal performance. Thus, the primary seal cannot, in 
general, be visualized as two perfectly flat and parallel surfaces. Some possible geometries 
are illustrated in Figure l-4[2]. The waviness (geometry a) and angular misalignment 
(geometry b) are most likely sources of hydrodynamic pressure build-up. Coning 
(geometry c) affects the hydrostatic pressure distribution and film stiffness. Externally 
imposed axial vibration (geometry d) can produce squeeze film damping. Parallel 
misalignment (also called radial eccentricity) and shaft whirl (geometry e & f) impart a 
radial velocity component to the fluid particles which can affect the leakage. Some of these 
seal geometries, particularly the angular misalignment (geometry b) introduces dynamic 
forcing function with frequency same as the shaft rotation rate on the flexible ring which 
would consequently exhibit oscillations in axial and angular modes. One particular interest 
to the seal designers is to whether the flexible ring would be able to dynamically track the 
rotor without metal to metal contact, for a given amount of rotor misalignment (commonly 
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called 'runout'). In the 'dynamic tracking' analysis, the fluid film is modeled as nonlinear 
springs and dampers. The damping comes from the 'squeeze film effects. The 
elastomeric O-ring also provides considerable amount of stiffness and damping both of 
which are frequency dependent. If O-ring starts slipping, Coulomb damping also comes 
into play. It has been found that higher the film stiffness is, the better the dynamic tracking 
ability of the seal ring becomes. 

The engineering of seals involves fluid mechanics, heat transfer, elasticity, 
thermodynamics - equilibrium and nonequilibrium, statistical mechanics, dynamics, 
chemistry and metallurgy, to name a few of the most frequent areas of concern. Usually 
each effect can be analyzed by itself, but then the integrated effects must be evaluated for a 
complete analysis of a sealing system. As indicated previously, seals are characterized by 
surfaces in relative motion separated by a very narrow gap. In order to ensure proper 
operation, very small differences in the dimensions of the seal part must be maintained 
while in operation. Deformations in geometry due to imposed thermal gradients, frictional 
heating, pressure and mechanical contact forces must be held to a minimum. In any case, 
net deformation must be no more than microvalues. Depending on the imposed conditions, 
seals operate basically in three different lubrication flow regimes shown in Figure l-5[2]. 
The 'non-contacting' seals usually operate with 'full film lubrication' whereas the 
'contacting' seals can operate either in the 'boundary (also called mixed friction ) 
lubrication' or 'dry sliding lubrication' regime depending on the excess magnitude of the 
closing force over the available opening force. 

If the sealed fluid is a gas (usually operating in a noncontacting mode), sometimes the 
'mean free path’ of the molecules may be of the same order or more than the nominal 
operating seal clearance (Knudsen No. > 1) in which case continuum fluid mechanics with 



no-slip boundary conditions is no longer valid and ’slip flow theory’ and 'statistical 
mechanics' are needed to describe the fluid flow. 

The selection of materials for seal rings is also a very important aspect of seal design. 
It requires extensive tribological testing to come up with a suitable material combinations 
for a certain specific application. In general, the seal ring materials should have good 
mechanical and thermal shock resistance, wear characteristics, corrosion resistance, self 
lubrication property, high modulus of elasticity, to name just the important ones. Carbon- 
graphite usually meets most of the requirements. It is quite frequently chosen in 
combination with some other compatible hard material like tungsten or silicon carbide. 
There are other combinations of seal materials used also, e.g., carbon-graphite vs. stainless 
steel, tungsten carbide vs. tungsten carbide, etc., depending on the operating conditions. 

With all the complexities and highly coupled effects that govern a seal behavior, it is 
not a secret that reliable and accurate design analysis for face seals does not exist. Any new 
seal design must be tested in the laboratory because a prediction of eventual performance is 
not possible on purely theoretical basis. Again there is a wide variation in performance of 
seals of "identical design;" a particular seal may fail after a few hours whereas another seal 
belonging to the same class can last for several years or so. It is because of this reason that 
seals are often termed the 'most unpredictable machine element' used in industry. 
However good a design and analysis tool is quite useful in evaluating one design against 
others. This procedure eliminates the need for building expensive prototypes and running 
time intensive laboratory tests for those designs which seem not so viable at the analysis 
phase [4], Also modeling and analysis give more insight into the complex mechanism of 
seal behavior. Hence there has been quite a bit of analytical and experimental work done in 
the face seal area over the last 25 years and efforts are constantly being made by the 
engineers and scientists to come up with better theoretical models for seal operation. Since 
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it is extremely difficult to perform a comprehensive seal analysis including all different 
coupled effects, most analytical work has been focused on one or two aspects. As of yet 
some of the individual effects (e.g., two phase flow modeling, high Knudsen number 
flow, mixed friction regime, wear model, nonequilibrium effects, to name a few) are not 
fully understood. 

1.1 Liquid Seals 

A number of investigators have analyzed face seals operating with incompressible 
liquid for different flow geometries. Etsion has extensively studied the angular 
misalignment effects on seal performance and stability. A misaligned face seal is shown 
schematically in Figure 1-6 [11]. He obtained a complete system of forces and moments 
acting on the flexible ring for different values of angular misalignment These can then be 
used in a seal dynamic tracking analysis. Etsion [5] observed that any angular 
misalignment produces a radial force on the flexible ring which in turn causes a radial 
eccentricity. When this eccentricity is large enough, the pumping of fluid may take place 
which will affect the leakage. The seal coning, however, tends to reduce the magnitude of 
the radial force [6]. When the pumping takes place in a direction opposite to the hydrostatic 
pressure drop, it is known as 'inward pumping.' This phenomenon was studied both 
analytically and experimentally by Findlay [7, 8]. Analysis also showed that a flat outside 
pressurized seal with angular misalignment has negative axial and angular stiffness [9]. 
Also the hydrodynamic forces creates a transverse moment which leads the tilting moment 
by 90 degrees [10] that can cause seal wobble. However, with coning the stiffnesses 
might change sign depending on the relation between the angle of tilt and angle of coning 
[11, 12]. For noncavitating flow, the effect of coning reduces the hydrodynamic 
transverse moment which would improve seal stability. The 'narrow seal approximation 
is usually made in seal analysis for simplification. With this approximation, the 
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circumferential pressure gradient and seal curvature can be neglected. Etsion [13] 
compared the accurate results from numerical solutions with the approximate results and 
found that over a radius ratio, r,/r 0 , greater than 0.8, an accuracy to within 1% can be 
obtained and hence in most cases this approximation is justified. 

For low pressure and/or high speed seals, lubricant cavitation is possible due to 
hydrodynamic effects. This has been experimentally observed. An interesting work on 
this subject has been published by Findlay [14]. The lubricant cavitation helps generating 
extra opening force because it prevents the generation of hydrodynamic pressure below the 
local vapor pressure of the liquid while not restricting the upper bound of the pressure. If 
cavitation did not occur, the components of hydrodynamic force would usually balance out 
and no net increase over the hydrostatic force would exist, which is not the case for low 
pressure seals. 

Sneck was one of the early investigators who made a very important contribution in 
the face seal analysis under incompressible flow. He published a series of papers [15] 
through [20] in '68 - '69 in which he addressed different aspects affecting seal 
performances, e.g., angular misalignment, radial eccentricity, tangential waviness, flow 
turbulence, centrifugal inertia and thermal effects. The centrifugal inertia term is included 
in the misalignment analysis in [15]. The centrifugal effects are shown to play a significant 
role in seal performance at higher speeds. For an outside pressurized seal, the regions of 
flow field may exist with a radially inward flow along the stationary surface and outward 
along the rotating one and under certain circumstances, there can be net zero leakage. The 
existence of such a region is a direct consequence of centrifugal inertial effects. This 
reverse flow phenomenon has been studied in detail in another paper [16]. The combined 
effects of misalignment and radial eccentricity is presented in [19]. The resulting leakage 
component can be outward (opposite to the direction of hydrostatic pressure drop) or 
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inward (in the same direction as the hydrostatic pressure drop) depending on the phase 
angle between the misalignment and radial eccentricity. Sneck also studied eccentricity 
combined with surface waviness [20]. Again the direction of leakage component is shown 
to be dependent on the phase angle. The once per revolution waviness is found to be the 
main contributor in the pumping effect. Turbulent flow is analyzed in [17]. The turbulent 
nature of the flow is described by an isotropic apparent viscosity model and a power law 
velocity profile. The misalignment and surface waviness are found to be somewhat less 
influential with turbulent flow than with laminar flow. 

The analysis of face seals is often based on the isothermal flow assumption within the 
seal clearance. The validity of this assumption is usually argued on the basis that seal faces 
are often good thermal conductors and hence will not permit large radial temperature 
variations. But even when the seal operates approximately isothermally, the temperature 
within the seal clearance need not necessarily be same as the cavity fluid temperature. An 
accurate prediction of seal performance requires an accurate evaluation of the fluid viscosity 
within the clearance space. A general thermal analysis procedure is presented in [18] to 
estimate the fluid operating temperature level inside the seal. No attempt has been made 
here to model the heat conduction through the seal rings. The upper and lower bound on 
the operating temperature are obtained by assuming adiabatic wall condition and zero 
thermal convection by the fluid, respectively. In a recent review paper by Khonsari [21], 
an extensive survey of literatures pertaining to thermal effects in slider and thrust bearings 
is presented with summary of important contributors of leading researchers and designers. 
Since thrust bearings and seals have some similarity, this paper is referred here. One very 
common assumption made by the seal analysts is to neglect the fluid temperature variations 
across the film. But viscosity variation across the lubricant film has been sometimes found 
to be responsible for generation of an appreciable load. In [21] many papers are cited 
which indicate the importance of transverse viscosity variation. King and Lauer [22] 
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presented an experimental method by infrared spectroscopy to verify the existence of the 
temperature gradients through the film. 

Pinkus and Lund [23] also considered the effects of centrifugal forces in high speed 
seals. They mentioned that at the upper limits of laminar conditions, centrifugal forces 
reduce the load capacity considerably and alter the pattern of the lubricant flow. Koga & 
Fujita [24] included both the radial and centrifugal inertia terms in their analysis of high 
pressure water pump seals. They obtained better correlation of the analytical predictions 
and experimental results when inertia effects are considered than their previous analysis 
neglecting these effects. 

As mentioned before, the total closing force is supported by hydrostatic and 
hydrodynamic fluid pressures and often by partial contact of seal faces (in contacting mode 
of operation only). For moderate to high pressure applications, the hydrostatic force 
component is predominant over the hydrodynamic component [25]. Since the film 
thickness is usually very small (of the order of a few microns), any local surface 
deformations due to the interfacial pressure and the angular twist of the seal rings under 
pressure strongly influence the hydrostatic load support and hence the seal performance 
[26]. For carbon rings with a relatively low modulus of elasticity, the distortions can easily 
be of the same order of magnitude as the nominal clearance of the seal. Thermal distortions 
can also occur due to both axial and radial temperature gradients in the seal rings caused by 
the frictional heat generated at the interface. Any radial taper in the direction of the flow 
changes the hydrostatic pressure distribution and the film stiffness. A diverging seals (in 
the leakage direction) exhibits a negative axial stiffness which may lead to seal collapse and 
high wear. 
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With the availability of digital computers and necessary softwares, the finite element 
(FE) analysis is commonly used for accurate predictions of pressure deformation and 
thermal distortion. A very important series of papers [27, 28, 29, to mention a few of 
them] have been published over the years by Metcalfe and his research staff at the Atomic 
Energy of Canada Ltd. (AECL) describing these analytical techniques. Analysis of seal- 
ring deflections due to applied pressure loadings, to thermal effects and to Coulomb friction 
between components is described in [28]. The deflection sensitivity of seal components is 
expressed as 'influence coefficients,’ evaluated with finite element analysis. He noted that 
Coulomb friction gives rise to undesirable performance hysterisis when operating 
conditions are changed. The correlation of experimental results and theoretical predictions 
is presented in [29], Salant [30] presented an analytical model of a generalized mechanical 
seal incorporating the fluid dynamics of the film and the mechanical and thermal distortions 
of various seal components. He utilized the concept of 'influence coefficients used in 
[28]. He found that the hydrodynamic forces due to waviness, roughness, misalignment 
and eccentricity produce insignificant opening force effects in comparison with the available 
closing force for high pressure seals. The hydrostatic pressure is responsible in carrying 
most of the applied load. Hence the face deformation, particularly 'coning' is the most 
likely controlling mechanism for load support for these kinds of seals. Based on this idea, 
Salant presented a novel design of an electronically controlled seal in [31]. A 
microcomputer based real time control system and electro-mechanical actuator dynamically 
adjust the seal coning and hence the film thickness, based on information received from the 
stator which monitors conditions of the film. This arrangement can greatly reduce face 
contact while limiting leakage by continuously optimizing film thickness. This would lead 
to a reduction in seal damage and wear and increase in seal life. 

Li [32] presented a finite difference heat conduction model for calculating the 
temperature distribution in the seal rings and resulting deformations. He considered one 
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dimensional toroidal deformation model in which the seal rings will preserve the geometry 
of the radial cross section after deformations. Doust and Parmar [33, 34] numerically 
analyzed axisymmetric distortions due to the pressure and thermal effects using 'boundary 
integral element' (BIE) method and correlated the results with the experimental 
measurements. They remarked that the BIE method is substantially more economical in 
terms of both computing time and storage than FE for the same level of accuracy. The 
main object of their test was to measure the fluid film geometry using capacitance type 
proximity probe, as a function of pressure. The sealant pressure and thermal effects 
essentially caused toroidal rotation of the faces for those seals used by them. Their rotation 
rate was found to be fairly insensitive to the interface pressure profile. They also observed 
hysterisis effect due to secondary seal friction. In a recent paper [35] by them, the effects 
of thermoelastic transients have been presented. The transient thermal distortion can be an 
order of magnitude greater than that at steady state. Transient response is worse for a 
shorter section than a longer one, although the time to reach steady state can be more than 
an hour for a long seal component. This is an interesting work since field surveys do 
suggest that transient operations can be more detrimental to the seal life than steady state 
running. 

Not as much work has been done in the 'mixed-friction lubrication' area for the 
contacting seals as in the 'full-film lubrication’ regime for the noncontacting mode of 
operation. The obstacle to further advancement in contacting seal technology is that 
relationship between controllable design parameters and performance parameters are not 
well understood. Lebeck has published a number of papers [36, 37, 38, 91] on 'mixed- 
friction' flow modeling and contacting seal analysis. He developed a model [36] which 
takes into account load sharing between mechanical and fluid hydrostatic pressure. The 
effect of wear is also modeled in order to predict how the radial profile alters and influences 
the hydrostatic pressure distribution with time. The experimental evaluation of the model is 



reported in [37]. In a contacting seal operation, an unstable phenomenon is observed by a 
number of researchers, including Kennedy and Grim [39], in which case a very slight 
amount of initial waviness on seal faces grows during seal operation. When this unstable 
condition, called ’thermoelastic instability' occurs in an operating face seal, the 
consequences - nonuniform wear, accentuated waviness, and high localized stresses and 
temperatures - can be very detrimental to seal performances. Kiryu et al. [40, 41] reported 
the generation of a "ringing" sound in a contacting water pump seal. They attributed this 
phenomenon to self-excited vibration due to ’stick-slip’ action, caused by transferring from 
fluid lubrication to dry sliding condition. Vibration mode in ringing sound generation is 
found to be mainly caused by the torsional and axial vibrations of the rotating shaft system. 


The previous investigations, mentioned so far, are mainly based on steady state 
analyses. However, the angular misalignment is inevitably present on the rotating ring 
which introduces a dynamic forcing function on the flexibly mounted stator. Hence the 
ability of the stator ring to track the rotor in a controlled manner is of great importance for 
safe seal operation and as the demand for higher operating speeds in rotating machinery 
increases, the importance of seal dynamics becomes more and more evident. Several 
researchers, namely Etsion, Green, Metcalfe and others, have addressed this issue 
analytically and experimentally in [42] through [55]. A review of face seal dynamics 
covering the literature until 1981 is presented in [50]. 

The flexibly mounted stator has basically three major degrees of freedom - axial and 
angular about any two orthogonal diameters. The twisting motion about axial direction is 
prevented by antirotation locks. If the radial stiffness of the O-ring secondary seal is low, 
which is usually not the case, then the stator can also move in the two perpendicular radial 
directions. The rotor transmits its angular motion to the stator via the thin fluid film 



separating the two seal rings. For a given forcing function, the response of the stator 
depends on its own inertia, the stiffness and damping of the fluid film and the elastomeric 
O-ring. The fluid film damping comes from the squeeze effects. In some cases the 
squeeze effects are an order of magnitude higher than the combined hydrodynamic and 
hydrostatic effects and hence play an important role in dynamic behavior of face seals. The 
stiffness and damping coefficients of the fluid film, both direct and cross-coupled, in the 
three major d.o.f. are calculated in [44, 51] based on small perturbation theory. It has been 
found that the narrower the seal, the less is the damping coefficients and at very small tilt, 
translation and rotational direct damping coefficients are an order of magnitude higher than 
the cross-coupled ones. The damping and stiffness characteristics of elastomeric O-ring are 
dependent on the amplitude and frequency of excitation and amount of squeeze. The 
experimental determination of the O-ring dynamic properties are presented in [54, 56, 57]. 
With a large rotor misalignment, the stator response is usually large and sometimes sliding 
and takes place at the O-ring interfaces and then the Coulomb friction becomes important. 

Dynamic analysis [45, 46, 48, 52] based on linearized small perturbation theory 
revealed three modes of operation: a stable mode in which a misaligned rotor is 

synchronously tracked by the flexibly mounted stator; a transition mode in which half- 
frequency wobble of the stator is superimposed on the previous synchronous tracking 
mode; and an unstable mode characterized by uncontrolled vibration of the stator, 
eventually causing failure. In the unstable mode, a seal will fail even with zero rotor 
runout. For low and moderate speeds, the stable mode seems to predominate. The stator 
tilt, however, differs from that of the rotor both in magnitude and direction. The difference 
and phase shift between two tilts result in relative angular misalignment between the rotor 
and stator. If this relative misalignment becomes too large, seal failure due to excessive 
leakage or even rubbing contact can occur even though the seal is dynamically stable. In 
[55], the complete nonlinear equations of motion of the stator are solved numerically. The 
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assembly tolerances in the form of initial stator misalignment and the dynamic properties of 
the elastomeric O-ring are accounted for in the analysis. Both stability threshold and 
steady- state response of the stator are investigated. In general, it was found that the critical 
shaft speed corresponding to stability threshold is quite high. Hence, the dynamic stability 
should not be a problem in the majority of noncontacting seals. A more practical problem 
is the steady-state dynamic response of the stator resulting from rotor runout and assembly 
tolerances. The results of the numerical analysis were compared with those of the previous 
small perturbation analysis that provides much simpler closed form analytical solution. 
Very good correlation was found between the two analyses for most cases of practical 
applications. 

Etsion and Burton [43] observed self-excited oscillations of seal ring in the form of 
precession and nutation. The wobble frequency was measured to be about 43% of the 
rotational frequency. Metcalfe [49] analyzed and tested a well-aligned face seal. He found 
that if the balance ratio is below a certain critical value the seal becomes hydrostatically 
unstable. If the elastomer stiffness in the tilting mode is insufficient to overcome this 
hydrostatic instability, the stator will exhibit wobble motion. The precession rate is 
theoretically found to be half the shaft speed if elastomer damping is insignificant (pure 
"whirl") and progressively slower as damping increases. Etsion presented an experimental 
observation of the dynamic behavior of face seals in [53]. The forced response of the 
stator due to the rotor runout was monitored by means of three proximity probes. It was 
found that both the stator tilt and its phase shift with respect to rotor tilt are time dependent 
and vary synchronously with the rotor rotation. The time variation is attributed to the 
presence of two components of stator tilt. One component is fixed in magnitude but tracks 
the rotor tilt. The other component is fixed both in magnitude and direction and is due to 
nonaxisymmetric effects in the flexible support of the stator. As a result, the relative 
misalignment between the stator and rotor was found to be time dependent. The 
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dynamically unstable seal behavior was also observed. At low supply pressure which 
means low film stiffness, result showed a sinusoidal perturbation at double the shaft 
frequency superimposed on the initial wobble due to angular misalignment. In the previous 
experiments, perturbation was always observed at half rather than double the shaft rotation 
rate. As the supply pressure was increased, this double frequency stable to unstable 
transition became a half-frequency transition instead. This higher frequency instability was 
not fully understood. 

1.2 Gas Seals 

The efforts on gas seal development started a little later than the liquid seals, and the 
work in this area are not so voluminous. In earlier times, the machinery, having gases as 
working fluids, like gas compressors, used, and some of them still use, liquid seals with 
an oil-buffered arrangement. The reason for this is that proper technology was not 
available to insure a non-contacting mode of operation, which is absolutely essential for gas 
seals because of the poor lubrication properties and high speed of operation. The oil is kept 
at a pressure a little higher than the sealed gas to ensure that only oil leakage could take 
place into the gas and seal would never run dry. In addition oil also leaks to atmosphere 
through another seal. Apart from the cost factor (about two orders of magnitude higher 
than the corresponding single phase seals), this design has some major disadvantages in 
terms of auxiliary equipment and space requirements. Also, since product contamination 
with just a small amount of buffer fluid may create enough problems, contacting liquid 
seals are typically used, which have inherently low and unpredictable life. Hence the need 
for a noncontacting seal development did arise for sealing gaseous fluids. 

Although the basic concepts are the same, the main difference between the liquid and 
gas seal analyses is that the governing equations describing the gas flow are nonlinear 
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because of compressibility effects and inclusion of flow inertia, which are sometimes too 
important to ignore. The flow is often turbulent and choking may occur at the outlet. Also 
under conditions of high velocity, the entrance loss effects cannot be neglected. The gas 
seals should also operate with a high film stiffness in order to have good dynamic tracking 
ability to prevent contact. As mentioned before, the pressure within the fluid film is 
generated hydrodynamically by the relative motion between uneven sealing surfaces and 
hydrostatically by frictional pressure drop through the seal. The hydrodynamic action 
ceases when the motion stops. There is no hydrodynamic pressure generation with parallel 
faces. To a limited extent all seals possess some hydrodynamic characteristics as a 
consequence of geometric imperfections and unplanned unevenness such as inherent or 
pressure induced circumferential waviness or micro-irregularities. These effects are usually 
quite small. The hydrostatic effects alone impart zero stiffness to a seal unless there is a 
radial coning in the flow direction. Because of these facts, some conscious efforts have 
been made to enhance the hydrodynamic action rather than rely on chance variation, by 
having planned uneven hydrodynamic patterns on the seal surfaces. Some of the 
commonly used patterns are spiral groove, Rayleigh-step pads, radial grooves, as shown in 
Figure 1-7 [58]. These are called 'hybrid' seals. 

The hydrodynamic pattern is followed by a seal dam which offers restrictions to the 
fluid flow and most of the pressure drop takes place there. Because of hydrodynamic 
action, there are some areas of higher pressure and other areas with lower pressure. Figure 
1-8 [58] shows the elevated pressure areas on the two seals. Figure 1-9 shows the 
components of 'F-h' curves for the hydrodynamic and the hydrostatic sections of a 
'hybrid' seal with Rayleigh-step pads, analytically obtained by Shapiro [59]. The two 
curves must be combined to get the net film characteristics for the seal under consideration. 
No angular misalignment effect is considered in this analysis. It is evident from this figure 
that the hydrodynamic action indeed imparts a very high film stiffness, particularly at small 
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clearances and prevents the seal faces from touching each other. The closing force is 
usually chosen so that the seal operates near the high stiffness region. Also it is seen that 
hydrostatic stiffness is almost zero and hence it does not contribute to the seal stability, 
although it may carry a major part of the closing force. Experiments performed by Ludwig 
[60] showed that seals with hydrodynamic pads outperformed the conventional seals used 
in small gas turbines. 

Some of the important research work on gas seals is documented in [59] through 
[72]. Cheng [61] analyzed a few different designs and found that the spiral groove design 
gives higher stiffness than Rayleigh-step pad one. This same conclusion was also drawn 
by Sedy [58]. He also brought out an interesting point. As mentioned before, most of the 
gas expansion takes place over the dam. The cooling effect associated with the expansion 
is sometimes several times more than the heating effect due to viscous dissipation. The net 
effect is the cooling of the gas near the dam and consequently a considerable amount of heat 
conduction takes place from the seal rings to the gas in the vicinity of the dam. The 
temperature gradient, thus set up in seal rings, tends to distort the seal face in a way to 
produce a divergent flow passage which has an unstable effect and sometimes causes seal 
contact at the outer diameter. Sedy suggested a wider dam design to overcome this 
problem because a wider dam would cause a higher heat generation which in effect tends to 
neutralize the cooling effect due to gas expansion. Zuk [64] presented a quasi one- 
dimensional analysis for the flow of gas through seals. This model includes fluid inertia 
and entrance losses, in addition to viscous friction which is accounted for by a friction 
factor. Subsonic and choked flow conditions have been predicted and analyzed. This 
model is valid for both laminar and turbulent flows. Hsing and Carraro [74] used an 
efficient algorithm based on fourth order Runge-Kutta with adaptive step size to solve the 
same governing differential equation. Shapiro [59] performed both steady state and 
dynamic analyses of a gas seal for jet engines. The seal dynamic response was found as a 
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function of rotor misalignment and secondary seal friction. He theoretically obtained 
superharmonic response, about four times the shaft RPM, which has been confirmed 
experimentally. This phenomenon has been attributed to nonlinear characteristics of the O- 
ring. The flexible ring is found to lose its tracking ability if the rotor runout or the friction 

force is too large. 

If the gas inside a seal is at sufficiently low pressure, the molecular mean free path 
can become comparable to the film thickness. The fluid subjected to this condition does not 
behave entirely as a continuum fluid but rather exhibits some characteristics of molecular 
chaos. One may also expect to encounter these effects in regions having very sharp 
gradients of fluid properties such that these properties change appreciably in the space of a 
few mean free paths, regardless of whether or not the absolute density of the gas flow is 
especially low. The dimensionless ration, A /h (Knudsen number), is a measure of the 
degree of rarefaction. When this ratio is large, the flow phenomena are mostly dictated by 
the molecular-surface interaction. This class of fluid flow is defined as free-molecular 
flow." For flows in which the value of Knudsen number is small, typically 
0 01 < / //i ~ 0. 1 , but not negligible as those in continuum mechanics, some departures 
from the usual continuum flow phenomena may be expected to occur. The layer of gas 
immediately adjacent to the solid surface no longer assumes the same kinematic condition 
as the solid boundaries but has a finite relative "slip velocity" and hence produces an 
apparent diminution in fluid viscosity. This is called the "slip-flow regime." 


The few research works covering these non-continuum effects in lubrication are 
documented in [73] through [78]. Hsing and Malanoski [74] found that if the lubricant is 
one of the gases having a large molecular mean free path, such as Helium, Neon or 
Hydrogen, the slip-flow phenomena could contribute substantial reduction in the 
performance of a thrust bearings, which is quite similar to face seals. Gans [75] derived a 
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slip flow lubrication equation for an arbitrary Knudsen number from kinetic theory. Fukui 
and Kaneko [76] developed a more accurate generalized lubrication equation based on 
linearized Boltzmann’s equation. The experimental results obtained by them [77] agreed 
well with their numerical results. Kubo et al. presented a finite element solution of the 
Boltzmann's equation in [78]. 

1.3 Two-phase Seals 

When liquid is sealed at temperature higher than its saturation temperature at the outlet 
pressure, it flashes inside the seal due to the pressure drop and/or the viscous heat 
dissipation. Typical examples of applications where such two-phase flow may be 
encountered are light hydrocarbons in petroleum refineries, hot water in boiler feed pumps 
and reactor coolant pumps, and cryogenic fluids like liquid oxygen and hydrogen 
(LOX/LH2) in rocket turbopumps. The two-phase seals generally exhibit more erratic 
behavior than their single phase counterparts. The seals also have more stringent 
requirements in their performances because of severity in applications. As, for example, 
light hydrocarbons are potentially flammable and explosive and hence certainly dangerous 
if allowed to leak. Since these hydrocarbons in gaseous phase are heavier than air, they 
usually form a thick dense cloud on the ground around the source. It constitutes a severe 
hazard [79]. In LOX/LH2 turbopumps, any seal failure due to excessive leakage can be, 
needless to say, extremely dangerous. Actually the face seals in the space shuttle 
turbopumps failed repeatedly on the test pads until they had been replaced with annular 
seals. Although annular seals are safer in operation, they allow very high leakage. At a 
later date, the face seals have been adopted successfully in the LOX/LH2 turbopumps for 
the Japanese H-l rocket [80]. Two-phase seal operation is also encountered in boiler feed 
pumps. It has been estimated that the boiler feed pump outages alone cost power 
companies several hundred million dollars each year in lost power revenues. It is believed 
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that a high percentage of these problems is related to seal failures [97]. The reactor coolant 
pump (RCP) seals can also experience change of phase of the sealed fluid during station 
blackout conditions and exhibit excessive leakage. Failure of these precision components 
may result in a small loss of coolant accident (LOCA) in nuclear reactors [82]. 

Because of the severity of application, the two-phase operation mechanisms must be 
better understood in order to come up with a suitable design. The research works done in 
this area are reported in [79] through [102], An interesting earlier paper on this subject was 
published by Orcutt [83]. He used a quartz runner to permit visual observation of the seal 
interface during operation. The experimental observations indicated the existence of a 
multiple phase film, characterized by two large scale regions. The first region adjacent to 
the seal cavity was occupied almost entirely by water. The second annular region extends 
from the atmospheric edge of the interface to a semi-stable boundary with the liquid-filled 
region. This region was occupied by a mixture of liquid and vapor. The boundary moved 
towards the edge adjacent to the seal cavity with the rise in liquid and seal surface 
temperatures. Unstable operation was encountered with visible leakage as the cavity fluid 
temperature was increased. More than a decade later, Harrison and Watkins [84] and 
Wallace [79] reported a similar unstable two-phase operation with light petroleum products 
at elevated temperature. Under the unstable operation, the fluid film periodically broke 
down and reformed with violent fluctuations in torque. Seals showed both audible (while 
in operation) and visible (when taken apart) signs of distress. Seal operation was, 
however, stable at lower temperature. Barnard and Weir [85] reported seals operating 
successfully with no visible leakage because of vaporization. The seal faces, they 
examined, all exhibited three concentric bands across their surfaces. Will [86] also 
observed the similar three banded appearance on successfully operating two-phase seals. 
No convincing causes are known. 
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In theoretical modeling of two-phase flow through seals, most of the work so far has 
been done by Hughes and his coworkers and are reported in [87] through [93]. Basically 
two different models have been presented in [87, 88, 89] for low and high leakage flows, 
respectively. The low leakage flow model is based on laminar flow and it considers heat 
conduction into the seal rings but neglects heat convection by the fluid and heat of 
vaporization. Boiling is assumed to be taking place at a discrete interface. The high 
leakage flow model is based on turbulent flow. This model disregards heat conduction 
through seal rings but takes into account convection, heat of vaporization, radial and 
centrifugal inertias. This model could predict continuous boiling over a finite region and 
also choking at the outlet under certain conditions. Beatty and Hughes [90] refined the 
turbulent flow model with better treatment of inlet losses. They obtained an anomalous 
'all-liquid choking' situation in which the flow is choked but remains liquid all the way up 
to the seal exit. Beyond the exit, the liquid flashes immediately into vapor. 

Lebeck presented a mixed-friction model with phase change in [91]. He modified the 
flow equation for roughness effects and considered the load support due to asperity 
contacts. Hughes and his coworkers assumed an idealized semi-infinite heat conduction 
model, whereas Lebeck used a more realistic seal geometry and boundary conditions and 
implemented a finite difference scheme to solve for seal ring face temperature distributions. 
Lebeck's model is valid for low leakage rates only. 

Beeler and Hughes [98] performed a dynamic analysis in the axial mode. They used 
the quasistatic 'F-h' curve obtained by using the adiabatic model to represent the fluid film. 
Squeeze film effects were ignored. With this limited dynamic model, they predicted self- 
sustained oscillations under certain conditions whereas failure due to metal-to-metal contact 
under other situations. Zuber and Dougherty [94] modeled the process of condensation 
and evaporation and derived a generalized lubrication equation. The two-phase region is 
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treated as a dispersed homogeneous flow with thermodynamic nonequilibnum between the 
vapor and the liquid. Although research efforts in this aspect are somewhat limited because 
of the difficulty involved, experiments showed that the effects of condensation and 
evaporation can become of primary importance in determining static and dynamic 
characteristics of saturated vapor bearings and seals. 

The 'F-h' curves for two-phase seal operation obtained by different investigators has 
a peculiar feature and is shown schematically in Figure 1-10. The positively sloped side of 
this curve implies a negative film stiffness whereas the negatively sloped side means 
positive film stiffness. For a given closing force, there can be two operating clearances 
with the smaller one giving rise to unstable operation and the larger one stable. 
Vaporization also seems to inhibit leakage. 

The two models for two-phase seal operation, developed by Beeler and Beatty 
[104,112], work reasonably well at the two extremes - very low leakage rates with 
convection neglected and very high leakage rates with conduction neglected. Both models 
break down as soon as the effect neglected in the respective model begins to become 
important. In actuality, most two-phase seal operations take place in the intermediate 
leakage range when both conduction and convection are important. A preliminary model is 
developed here to bridge the gap between the two previous models. This model, known as 
the 'Film Coefficient Model,’ is valid over the entire laminar flow regime unlike the earlier 
model developed by Beeler which only worked at the very low leakage rate end. The new 
model considers both conduction and convection and allows continuous boiling over an 
extended region whereas the earlier model which neglects convection always forces a 
discrete boiling interface and exhibits numerical instability as soon as leakage rate starts 
becoming a little higher. With the inclusion of turbulence and radial inertia effects, the 
applicability of the 'Film Coefficient model’ can be extended to high leakage rate end with 
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the ability to predict choking. Hence this model has the potential for describing the seal 
behavior over the entire range of possible leakage rates - low to high. 


Another simplified and semi-analytical model, known as 'Isothermal Model, has 
also been developed for low leakage rates. This is based on the model developed by 
Beeler. The assumptions of isothermal condition along the seal interface and ideal gas 
behavior of the vapor permit closed form solutions which may be used for preliminary 
design and analysis. However, to obtain more accurate and realistic description, the Film 
Coefficient Model' may be used. 

Under certain two-phase operation, seals seem to exhibit self-sustained oscillations 
even when all the applied conditions remain quite steady. These have been observed as 
axisymmetric fluctuations in the film thickness accompanied with periodic interface 
temperature variations. 

1.4 Two-Phase Seals - How They Work 

We continue the general background on two-phase seals, laminar (low leakage) and 
turbulent (high leakage). Details of the equations and computational techniques will be 
presented in Chapter 2 and in more detail in the appendices, except for the simplified model 
which is discussed in detail in Chapter 3. 

Again, many phenomena associated with all liquid on all gas seals have been 
discussed in considerable detail and do provide insight into their behavior. However, 
many effects, such as popping, chattering, and some failure modes are associated with 
two-phase effects in that the behavior changes in response to temperature, subcooling of 
the sealed liquid and generally whether boiling occurs. 
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Liquid and gas face seals have, generally, neutral axial stiffness but can be unstable to 
wobble. In general, coning (convergence in the flow direction) tends to give positive axial 
stiffness to liquid seals and many are designed this way. 


However, two-phase seals can have a negative axial stiffness and be inherently 
unstable to axial disturbances over certain ranges of parameters. Coning may serve to 
mitigate the negative stiffness, but not always. In fact, a stably operating seal may become 
unstable by changing the operating conditions, particularly by increasing the temperature of 
the sealed liquid (i.e., decreasing the subcooling and approaching saturation conditions). 
These observations apply both to laminar and turbulent seals. 

In order to understand the characteristics of a two-phase face seal let us consider the 
flow through the seal. Figure 1-11 shows the trajectory in a T-s plane for flow through a 
seal. The actual distance from f to g on the seal face can vary from a negligible distance to 
the entire seal face. For most low leakage seals the points will be close together and boiling 
takes place almost at a discrete radius. For high leakage turbulent seals the boiling may 
occur over the entire seal face. Further, the closer to isothermal operation, the shorter the 
region over which boiling occurs. Clearly if there is no temperature change then the boiling 
must occur at a discrete interface in order to satisfy both momentum and the Clapeyron 
relation. 

Now consider the pressure drop through the seal. Figure 1-12. If the seal is all 
liquid, the pressure is nearly linear, and if all gas the pressure is nearly quadratic. 
However, if boiling occurs, then for a given Film thickness (seal face separation), the 
leakage rate is reduced and the pressure is higher than for an all liquid on all gas seal. A 
plot of the total opening force produced by this pressure vs. the seal face separation then 
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produces the axial stiffness curve Figure 1-10, which is shown again in Figure 1-13 with 
the liquid and gas asymptotes. A positive slope here represents a negative axial stiffness 
(i.e., instability), and a negative slope represents a positive axial stiffness (stability). 
Figure 1.14 shows a set of actual curves generated from the simplified theory which will be 

discussed in Section 3. 

These axial stiffness curves are the key to the distinctive two-phase operation. 
Operation on right hand side of the curve in Figure 1-14 is stable to an axial disturbance but 
the left side is unstable. Now, the seal may be balanced for any level by changing the 
balance ratio. An arbitrary line is shown in Figure 1-14 as the balance point. At each 
speed there are two equilibrium points. The one with the larger film thickness is stable, the 
one with the smaller film thickness unstable. Depending on the location of the balance line 
the seal may be unstable for an all liquid seal, all gas, or both. For instance, if the balance 
were established at 1000 N, the seal would open if it were all gas and collapse if all liquid. 
At the balance shown, about 1250 N, the seal would collapse if the seal were either all gas 
or all liquid and relies on two-phase operation for stability. 

However, the situation is more complex. The behavior depends critically on the 
subcooling of the sealed liquid. As the sealed fluid nears saturation conditions the stiffness 
curve tends to become entirely positive in slope and the seal is totally unstable. Hence, a 
seal balanced properly at one level of operation may become unstable if the temperature is 
changed sufficiently. These considerations are critical in such situations as nuclear power 
plant "black-outs." Consider Figure 1-15. The saturation temperature is 453 K. At about 
440 K the curve shows a monotonic positive slope (unstable) and the seal would tend to 
pop open. Interestingly, if the seal temperature were very close to saturation the opening 
load tends to correspond to all gas and the seal would collapse. During a transient, opening 
would occur first with possible catastrophic consequences. 
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As opening continues, what happens as the leakage flow increases and become 
turbulent? We can answer that with the turbulent flow model which is discussed in detail m 
Chapters 4 and 5. Figures 1-16 and 1 1-7 show turbulent curves. Generally, the turbulent 
curves show the same trend and we conclude that an instability is exascerbated as the seal 
continues to open or collapse. 

To summarize our findings, 


1 . Codes developed for steady operation description, stiffness calculations, and 
stability analysis for annular and face seals - laminar or turbulent. 

2. Axial disturbances may create instabilities and possible self-sustained 
oscillations. 

3. All behavior critically dependent on heat transfer effects and viscous 
dissipation. 

Before we review the general equations we use in the various models, a word about 
some of the necessary considerations is in order. We must consider the equations of 
momentum and energy, the thermodynamic equation of state, the Clapeyron equation, 
viscous dissipation in the fluid, centrifugal inertia in the fluid, and the heat transfer mto the 
seal faces. The consideration of heat transfer is crucial to the behavior of the seal. In 
Figure 1-18 we show a schematic drawing of how the heat can flow into or out of the fluid 

through the seal faces. 
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Figure 1-1: Classification Diagram - Industrial Sealing Devices 



Figure 1-2: Schematic Diagram of a Face Seal 
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Figure 1-3: Force Diagram on a Stator 
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Figure 1-4: Possible Primary-Seal Geometries 
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Figure 1-5: Seal Lubrication Models 
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Figure 1-7: Hydrodynamic Patterns on Face Seals 



Figure 1-8: Hydrodynamic Pressure Generation 
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Figure 1-9: Hybrid Gas Seal F-h Curves - Hydrodynamic and Hydrostatic 

Components 
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Figure 1-10: Seal F-h Curves for Two-phase Operation 
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Figure 1-13: STIFFNESS CURVES 
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Figure 1-15 
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Figure 1-16: 


Stiffness Curves Showing 
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CHAPTER 2 

GOVERNING EQUATIONS FOR LAMINAR SEALS 
2.1 Introduction 

The governing equations of fluid motion and heat transfer through the seal rings are 
developed here. Figure 2-1 shows the radial flow passage inside a seal under two-phase 
operation. Although this figure indicates a parallel face, the formulations are done with a 
more general coned face seal with film thickness, h, varying along the radial direction, r. 
The moving and the fixed plates represent the rotor and the stator, respectively. For this 
particular case, the heat generation is high enough and the flow is low enough to cause 
complete boiling inside the seal. There may be situations in which the sealed liquid does 
not boil completely and leaves the seal as a two-phase mixture. 

Heat is generated in the lubricating film due to viscous dissipation. A part of this heat 
is conducted into the seal rings and the rest convected downstream by the fluid. The 
proportion of each mode of heat removal at any radius depends mainly on the leakage rate, 
the material conductivities and whether or not phase change takes place. Usually under 
two-phase operation, conduction and convection both play important roles depending on 
the region under consideration - liquid, two-phase, or gas. 

The following assumptions are made in this analysis: 

1 For the heat conduction calculation, the seal geometry is to be modeled with the 
appropriate thermal boundary conditions. As the radial width of the interface is 
very small and backup materials on either side of the interface are large in 
volume, the seal rings are assumed to be semi-infinite solids of homogeneous 
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composition with conductivity same as the average value of the two materials. 
The idealized seal model is shown in Figure 2-2. The temperature very far from 
the interface, Too , corresponds to the bulk temperature of the liquid being 
sealed. These assumptions make the model independent of the actual 
geometrical details, and simplifies the numerical calculations. It is expected 
some general conclusions may be drawn about seal behavior using this model. 
For actual design purposes, however, the seal geometry and all the imposed 
thermal boundary conditions should be given proper consideration, and the 
methodology developed here can still be useful. 

2. The film geometry is assumed to be axisymmetric, with no tilting or angular 
misalignment. 

3. Only quasistatic, laminar flow is considered in the analysis. However, the 
governing equations for squeeze film effects, developed by Beeler [104] and 
turbulent, inertia dominated flow, developed by Beatty [90] are briefly 
discussed here for the sake of completeness. 

4 . The thermodynamic and transport properties are uniform across the film. 

5 . Density of liquid is constant throughout. This may be suitable for water but not 
for many hydrocarbons. 

6 . Radial inertia is neglected in laminar flow on the ground of low leakage rates. 
Although this is very small for low leakage seals, the speed of sound in a two- 
phase mixture, however, can be far less than that in either pure liquid or 
gaseous phase under the same conditions [105]. This usually happens around 
50% quality range. Hence 'choking' may occur even at apparently low leakage 
rates, radial inertia and inlet losses included in the development of the turbulent 
equations. 

7. Radial conduction of heat in the fluid film is neglected because of low 
conductivity of the fluid. 

8 . Centrifugal inertia effects are included. 
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2.2 Governing Equations of Fluid Motion 
2.2.1 Steady Laminar flow 


Momentum and Continuity 


Referring to Figure 2-1, the 0 momentum equations for the fluid are as follows. 


% + m ?JL = .P 
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( 2 - 1 ) 
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( 2 - 2 ) 


dz 


The term on the right-hand side of the Equation (2-1) represents the centrifugal inertia 
effect. 


Integrating (2-2) across the film twice, subject to the boundary conditions the w — 0 
at z = 0 and w = car at z = h. 


The Equation (2-3) is an example of the classical velocity driven Couette flow between two 
parallel plates. 


Substituting (2-3) into (2-1) and integrating with the boundary conditions u = 0 at z 
= 0 and h. 
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As seen in Equation (2-4), the radial velocity has two parts. The first term on the right 
hand side is the typical parabolic profile of Poiseuille flow. This is driven by the pressure 
gradient. The second term comes from the contribution of the centrifugal inertia. It has a 
stronger effect near the rotating face and distorts the Poiseuille profile which would 
otherwise be symmetric about the center of the film. 

The continuity equation in integrated form is 

r h 

m = p J 2nrudz (2-5) 

J 0 


where m is constant along the radius for steady-state seal operation. Substituting (2-4) in 
(2-5), there results. 
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The second term on the right-hand side of Equation (2-6) is due to the centrifugal inertia 
effects. The Equation (2-6) is valid for liquid, two-phase, or gas. 

Energy Equation 

The derivation of the integrated energy equation is shown in Appendix A. This equation 
can be written as follows: 
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The left hand side term represents the convection. The first term on the right hand side is 
the heat generation due to viscous dissipation. The second and third terms are the 
correction terms due to the retention of the centrifugal inertia in the equation of motion, and 
q is the heat conduction flux into the seal plates from the fluid film. 

Two-phase Region Model 

The two-phase fluid is assumed to be a homogeneous mixture of saturated liquid and 
saturated vapor. The specific volume, enthalpy and thermal conductivity are obtained by 
linearly weighting the corresponding saturated liquid and saturated vapor properties as 
follows: 


v = v / + ' lv / g 

(2-8) 

l=i f + u fs 

(2-9) 

k = kf+ Xkj g 

(2-10) 


The weighting parameter, X , is the quality or the local mass fraction of vapor. 

Often the details of the two-phase flow pattern are not known and an idealized 
theological model cannot be defined. Faced with the necessity of choosing some 
expression for viscosity, many workers have chosen averages which fit the limiting cases 
in which only either phase is present [105]. Two of such models are presented below: 

Model 1: Weighted according to the mass fraction, 

n=n f + Xn fg (2-11) 

Model 2: Weighted according to the volume fraction, 

ti =(l-X)^Hf+ ( 2 - 12 ) 
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which leads to the following expression for kinematic viscosity: 


With the Model 1, a two-phase region with very low quality will behave like a liquid region 
so far as the heat generation is concerned but the pressure drop will be rapid similar to a gas 
region, whereas with the Model 2, the heat generation will be sharply reduced over the 
same region but the pressure drop will be low as in liquid region. 

Specification of Fluid Properties 

Thermodynamic data and fluid transport property data is required to solve specific 
problems. To this end, saturation thermodynamic data tables were drawn from Reynolds’ 
book [25] and fit to fourth-order splines. Accurate interpolation between table entries could 
then be done. First derivatives of the saturation thermodynamic properties with respect to 
pressure are also required in some analyses and can be estimated with good accuracy from 
the spline curve equations. Saturated liquid and vapor viscosities for a number of 
substances of interest were found in a variety of sources [26, 19, 29]. Viscosities of 
superheated vapors were given in tabular form by the same sources and bilinear 
interpolation between those entries was done as necessary. 

Because the state of subcooled liquids and superheated vapors depends on a 
combination of temperature and pressure, calculation of fluid properties in those 
thermodynamic states can be unwieldy. However, these calculations can be simplified by 
some reasonably accurate approximations which are worthy of discussion: 

• The viscosity of a subcooled liquid was taken to be the viscosity of saturated 
liquid at the same temperature. In all cases under study, the given pressure levels 
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were only moderately above the saturation pressures. Since the dependence of 
viscosity on pressure is rather weak in all the substances considered, this 
approximation is judged to be adequate. 

• The specific volume for a subcooled liquid was likewise taken to have the value 
of saturated liquid at the same temperature. Liquids are, in general, highly 
incompressible so that this approximation is adequate as well. 

. The specific internal energy of a subcooled liquid differs little from the internal 
energy of a saturated liquid at the same temperature. As discussed in the previous 
item, the specific volume also changes little between these states, allowing the 
specific enthalpy of a subcooled liquid to be estimated with good accuracy. The 
subcooled liquid enthalpy is approximately equal to the saturated liquid enthalpy 
taken at the subcooled temperature plus the saturated liquid specific volume 
multiplied by the difference of the pressures of the subcooled and saturated states. 

• The need to calculate properties of superheated vapors is very seldom encountered 
in typical turbulent seal analysis. The means to make such calculations are 
included solely for completeness. Superheated vapor was considered to be an 
ideal gas with constant specific heats, but variations in gas viscosity with 
temperature and pressure were taken into account For thermodynamic states near 
the saturation dome, this approximation is judged to be rather crude. 

2.2.2 Transient Laminar flow & Squeeze-Film Effects 


Squeeze-film effects refer to changes in the fluid properties due to relative axial motion of 
the seal rings. As one ring approaches the other, the lubricating film is squeezed. Usually 
the pressure inside the seal is higher than the steady state solution with the same film 
thickness. The mass flow rate will in general vary along the radius. If the film thickness 
decreases quickly enough the pressure in part of the fluid film may rise above the sealed 
fluid pressure. If this happens, fluid will be forced out of the seal at the inlet as well as the 
exit [104]. On the other hand, when the seal rings move apart and the film thickness 
increases, the pressure in the fluid film will fall below that predicted by the steady state 
solution. Squeeze-film effects will always produce damping because the change m opening 
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load resulting from an axial velocity of the primary ring always acts in the direction 
opposite to that velocity. 


The fluid governing equations for single-phase flow, developed by Beeler [104], are 
given below. 


Momentum and Continuity - Single Phase 

Since radial and circumferential inertias are neglected, the r and z momentum equations 
will be the same as in the steady-state case. However, the pressure and the mass flow rate 
both vary along the radial direction as well as in time. Rewriting the Equations (2-3) and 
(2-4) for the circumferential and radial velocities, 
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* ph 
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(2-14) 


Following the derivation of the steady state flow, the Equation (2-14) is multiplied with the 
circumference and the density, and then integrated across the fluid film to get an expression 
for local leakage rate, m. 
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(2-15) 


Note that m may now vary with radial position. Considering the annular control volume in 

Figure 2-3, as the film thickness changes with time, the rate at which the mass within the 
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control volume changes must be equal to the difference between the leakage rates on each 
side of the control volume. This may be written mathematically as follows: 


m r - m r + sr = (2nrph8r) 


(2-16) 


Now dividing through by the circumference and the width of the control volume, 8 r , and 
then letting 5 r approach zero, there obtains, 

d(ph)_ _ 1 lim m r+fr' m r 

2n r 8r 

(2-17) 

1 dm 
2 nr ~^ r 


Finally, using the Equation (2-15), the following form of the squeeze-film equation of 
motion for single phase is obtained, 


pV + h 


dp _ J_ d_ 
dt ~ I2r dr 





(2-18) 


where V is the relative axial velocity of the two rings. 


Energy - Single Phase 


The energy equation is given below, with the terms neglected consistent with the earlier 
assumptions. 
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Equation (2-19) is integrated across the film and it can be written as follows: 


di dp po)r 


2 

ph 

2 4 
pr o) 

P 

112 




_/_ di dp 
12 Ik Hr 


~ Q 


( 2 - 20 ) 


Two-Phase Region 

Under dynamic situations, the two-phase region model is more complex than that under 
steady conditions. The appropriate field equations are developed by Zuber and Dougherty 
[94], The two phases may not be in thermal equilibrium and also the process of 
condensation and evaporation are to be modeled. Formation of condensate acts as a vapor 
sink, whereas the evaporation acts as a vapor source at a point. The effects of the vapor 
sink and/or the source term can be of primary importance in determining dynamic 
characteristics of two-phase seals [94]. Apart from the continuity equation for the mixture, 
which is same as the Equation (2-17) for the single phase, the continuity equation for the 
liquid is to be considered also, which is given below, 

+ V : ((1 - A) p V) = T f (2-21) 

where the source term Tf is the mass rate of liquid formation per unit volume. 


The following constitutive equation of condensation or of evaporation is also 
required, 


T f = T f (p) (2-22) 

For condensing flow, the liquid source depends on the rates of droplet nucleation and 
vapor condensation on these droplets. For evaporating flows, Tf depends on the droplet 
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number density, size distribution and the rate of evaporation. In either case, the 
constitutive equation of condensation and/or evaporation is a strong function of vapor 
saturation pressure. At present, the constitutive equation is not known for some flow 
regimes and only in a rudimentary manner for others. 


The amount of complexity involved in the study of two-phase seal transients has been 
presented here briefly. No attempt has been made to derive the full set of governing 
equations for these cases. 


2.3 Heat Conduction Through Seal Rings 
2.3.1 Steady-State Heat Conduction Model 


The heat transfer into the seal rings must be considered in order to determine the seal face 
(also referred to as 'wall') temperature distribution. In Figure 2-4, the face of one seal ring 
(both assumed to be a semi-infinite solid of thermal conductivity k) is shown. The 
temperature at position r may be expressed as 


T(r) = 



Q(r')dA 

______ p— — — 

4nk r-r' 


(2-23) 


where k is the average conductivity of the two seal ring materials and a is the annular seal 
interface over which heat is generated. 


Since the fluid film thickness is very small compared to the other relevant dimensions 
in the problem, it may be considered to be a thin heat "source" which releases or absorbs 
heat at a given rate, q , per unit area. The fluid film is modeled as an annular disk of 
infinitesimal thickness laying in a plane and divided into n coplanar annular elements of 
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equal width. This disk is imbedded in an infinite solid of thermal conductivity k which 
represents the seal ring and seal. Each element corresponds to one radial finite difference 
point which lies halfway between the element's inner and outer radii. With these 
assumptions, Equation (2-23) can be written in matrix form as follows: 


{T}= [<:]{<?} + 


(2-24) 


{T ) and {q ) are the temperature and heat generation vectors, respectively. is a 

column vector of constant elements. In this derivation, the heat dissipation q t is assumed 
to be constant over the element i . [C ] is called the 'steady-state influence coefficient 
matrix.’ The detail derivation of [C ] is done by Beeler [104] and the final expressions of 
its elements are shown in Appendix B. However, the elements of [C] can be obtained for 
finite seal geometries with proper boundary conditions by using either finite difference or 

finite element methods. 


2.3.2 Transient Heat Conduction Model 

Under certain circumstances, when the fluid film is not stable and exhibits variation in 
thickness over time, the seal interface temperature also becomes time-dependent to 
determine which the following semi-infinite transient heat conduction model is developed. 

The governing differential equation for transient heat conduction in three dimension is 


dt pC 


(2-25) 


where h is the heat generation rate per unit volume. 

The initial and infinity conditions are T(r,0) =f(r) and, as|r| ~> 00 , T --> 0. 
The Green's function for the Equation (2.25) is given by. 
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(2-26) 


G (r,t £,t) = 
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The Equation (2-26) gives the influence coefficient at a radius r at time t due to a unit 
source at £ at time 1 . 


The temperature at r at time t is given by, 


T <r,f) = 



G(r,t ;£,0 )pC /(£) <£ + 


f 0 t f~G(r,t;l,0)h(l,x)dldr + T, 


(2-27) 


The first term in the right hand side of the Equation (2-27) represents the contribution from 
the initial conditions, the second term from the heat generation and the third term comes 
from the steady background temperature. For a seal, the heat generation is taking place 
over the annular interface only. Hence while the spatial integration is to be done over the 
entire three-dimensional region in the first term, the same integration is only necessary over 
the finite interface region in the second term. One way of avoiding having to calculate the 
initial condition integral is to perform the time integration over the entire time interval (0, t) 
to find the temperature at the end of each discrete time step t n because at time t = 0, the 
initial temperature distribution f (<£ ) is zero everywhere. 


The fluid film is again modeled as an annular disk of infinitesimal thickness lying in a 
plane and divided into n coplanar annular elements of equal width, as is done in the steady 
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state case. At any time t , the entire time interval (0, t ) is discretized into N intervals of 
varying size. The disadvantage with this semi-infinite model, as mentioned above, is 
finding temperature distribution at t n , cannot be conveniently posed as an initial-value 
problem over the n th time interval with the known temperature distributions at the end of 
the previous interval tn-i as the prescribed initial condition. Instead the first term in 
Equation (2-27) is set to zero and time integration in the second term is performed over the 
entire time interval (0, t n ) which has been discretized into n sub-intervals. It is assumed 
that the heat generation vector {q } and temperature vector {T } are constant over each of 
these sub-intervals. With the chosen spatial and temporal discretization, the Equation (2- 
27) gives the following matrix equation with summation over all the time sub-intervals, 

{r(fj} = X[C*](<7*} + (r~) (2-28) 


where { T (t n )) is the temperature distribution vector at time t„ and [C* ] is the 'time- 
dependent influence coefficient matrix,' the detail derivation of which is given in Appendix 
C. 
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2.4 Boundary Conditions 

The seal inlet and outlet pressure p 0 and pi , and the bulk fluid temperature T„ are the 
given boundary conditions. The inlet pressure, p 0 , is the same as the sealed fluid pressure 
for low leakage rates. However, the inlet losses are to be accounted for high leakage rates. 
The pressure pi is same as the back pressure at the seal outlet (p „) except where choking 
occurs. The material temperature further away from the seal interface, Too , is assumed to 
be the same as the bulk sealed fluid temperature. Although a part of the heat generated at 
the seal interface is conducted into the seal rings and the surrounding materials, the bulk 
fluid temperature will not change much because the heat will be ultimately carried out of the 
machinery with the fluid flow. 
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Figure 2-1: Two-phase Flow Through Face Seal 
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Figure 2-2: Idealized Seal Model 
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CHAPTER 3 

A SIMPLIFIED MODEL FOR LAMINAR 
TWO-PHASE FACE SEAL DESIGN 

3.1 Introduction 

In this chapter we present a simplified model suitable for designing either parallel 
and coned face seals. During the early stages of designing a seal, it is usually necessary to 
perform extensive parametric studies before viable designs can be identified. To a seal 
designer, a computationally efficient model for quick and simple analyses is a very valuable 
tool. Fortunately, most face seals are designed to operate at very small film thicknesses 
and have very low leakage rates, additional simplifications to our basic laminar seal model 
presented in Chapter 2 can be made to greatly simplify the analysis. 

When the seal film thickness is very small, the fluid mass flow (leakage) rate is low 
enough that the heat carried away in the fluid (convection) is very small compared to the 
heat generated due to viscous dissipation. Also, because of the high heat generation rate, 
phase change of the fluid from liquid to vapor occurs over a very short distance. It is 
therefore reasonable to assume that, in a low leakage seal analysis, the errors introduced by 
neglecting the effects of convection and continuous boiling will not be unacceptable. In our 
simplified model for low leakage seals, the most important simplifications made are 
neglecting convection and assuming discrete boiling. A computer code suitable for PC 
operation based on this model has been prepared and is documented in Appendix D. 

3.2 Physical Model 

We will be concerned primarily with outside pressurized seals although the analysis 
is equally valid for inside seals. Figure 3.1 shows the geometry of a face seal. For an 
outside face seal, the leakage path is radially inwards. High pressure fluid enters the seal 
from the outer radius and leaves at the inner radius. This configuration offers the 
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advantage of a lower leakage rate with the help of the centrifugal effect. For an all liquid 
outside seal, net leakage can be reduced to zero if the seal is operating above a certain 
critical speed. However, in most operations, this critical speed is much higher than the 
operating speed and the inertial effects are of minor importance. 

In this simplified model, a number of physical assumptions in additional to those in 
Chapter 2 were made and a discussion of them is given below. 

1 . The seal is axisymmetric and the seal surfaces are perfectly smooth. There is no 
surface contact between the seal plates. For face seals operating at extremely small 
film thicknesses, hydrodynamic effects caused by surface roughness and surface 
contacts can contribute significantly to the net opening force. Since our model does 
not take these effects into account, the opening force estimate will err to the low 
side when such effects are indeed significant. However, coning is included in the 
model and wear at the outer radius may be approximated by an equivalent coning. 

2. The density and viscosity of liquid are assumed to be the same as saturated liquid at 
the same temperature and is independent of pressure. Viscosity of vapor is also 
assumed to be equal to that of saturated vapor at the same temperature. This 
assumption is good for water but is not very suitable for many hydrocarbons, 
where densities and viscosities are sensitive to pressure variation. Since no 
provision is made to handle fractional boiling, the equations developed are not 
applicable to multi -component fluids. 

3 . Centrifugal inertia effects in the vapor region are small compared with those of the 
liquid region and are neglected. 

4. The temperature of the fluid is assumed to be uniform across the film and is the 
same as the temperature of the surfaces of the seal plates. In reality, a temperature 
gradient must exist in the fluid film for heat transfer to take place between the fluid 
and the seal plates. However, the film thickness is so small that a very small 
temperature variation is sufficient to produce the required gradient 

5 The fluid flow is assumed to be quasi-isothermal at the boiling location temperature. 
When boiling is discrete and convection is negligible, determining the exact 
temperature profile of the fluid along the leakage path is not very critical to the 
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analysis. For low leakage seals, the variation in the fluid temperature along the seal 
faces is often small. The only relevant seal temperature is the temperature at the 
boiling location, where the boiling interface pressure is equal to the saturation 
pressure corresponding to this temperature and must be determined by looking up 
the saturation table. With the exception of saturation pressure, all the relevant 
physical properties, namely the fluid viscosity and density, are insensitive to small 
temperature changes. Under the assumption of negligible convection, viscosity is 
the only fluid property needed to calculate the boiling location temperature. The 
assumption of an isothermal flow, which in this case means constant viscosity of 
the fluid, will not introduce any significant error. Since the boiling interface 
pressure is directly evaluated using the boiling location temperature, the error 
introduced by the isothermal assumption is also minimal. After the boiling interface 
pressure is found, the leakage rate and the pressure profile of the seal are evaluated. 
The isothermal assumption, once again, will not introduce any significant error to 
the results because the only physical properties needed are viscosity and density. 

It should be noted that, the physical meaning of assuming a quasi-isothermal flow 
is that the temperature variation along the fluid leakage path is assumed not to be 
large enough to significantly affect the viscosity and density of the fluid. Since it is 
shown that no unacceptable error will be introduced, this assumption is justified. 

6. Boiling is assumed to occur at a discrete location. If the flow is indeed isothermal, 
then the boiling (change from liquid to vapor) must occur at a discrete interface. 
This is true because in any two-phase region the temperature and pressure are 
related by a unique relationship (Clasius-Clapeyron equation). To satisfy 
momentum balance, the pressure drop must be monotonic through the seal. Hence 
there can only be one position where the saturation condition is satisfied if the 
temperature is uniform. 

In reality, there must be some heat transfer "to" the fluid to effect boiling. For low 
leakage seals this amount of heat, and indeed the convection in general, is negligible 
compared to the total heat generated by viscous dissipation, which is essentially all 
conducted into the seal plates. 

7 The vapor is assumed to behave as an ideal gas. This assumption greatly simplifies 
the analysis and unless a more precise heat transfer analysis was used it would 
seem unwarranted to incorporate vapor properties of real fluids. Exact 
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thermodynamic data were used by Hughes and Chao [88] and quantitatively it was 
shown to make little difference. 

8. Since we assumed quasi-isothermal flow for the low leakage seal with discrete 
boiling, the model of a two-phase mixture becomes moot. We only consider pure 
liquid and vapor. However, it is possible that a separated flow may exist with 
vapor accumulating along the rotating face and liquid along the stationary face 
because centrifugal inertia would tend to sweep the liquid off the rotor. 

3.3 Mathematical Analysis 

For the liquid region, using the assumptions of constant density and viscosity, 
Equation (2.6) can be integrated to obtain a closed form solution of the pressure 
distribution. The boundary conditions applied are p = p 0 at r = r Q and p - p b at r = r b . 


For a parallel faced seal, the expressions for the mass flow rate, m, and the 


pressure, p, are, 
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For a seal with constant coning slope, (5, the expressions are. 
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From Equations (3.2) and (3.4), we can see that it is possible for the mass flow rate 
of an all-liquid seal to go to zero when the centrifugal inertia force balances the pressure 
force. The speed at which the mass flow rate is zero is called the critical speed and is given 
by the following expression. 

V 20 (p 0 -p) 

W'-rb ( 3 - 5 ) 

When the seal speed is at the critical speed, a liquid front will be formed at the exit. 
In theory, if the seal is running at above critical speed, the direction of the mass flow will 


be reversed (i.e. from the inner radius which has a lower pressure to the outer radius which 
has a higher pressure). However, it is not possible for most seals to have reversed flows 
because there is seldom any liquid available at the exit for being "pumped back into the 
reservoir. If the seal speed goes higher than the critical speed, the liquid front will move 
towards the inlet and significantly reduce the seal opening force. For a detail discussion on 
the effects of centrifugal inertia, the reader is referred to Basu [92]. 


For the vapor region we assumed ideal gas behavior, negligible inertia effects and 
vapor viscosity equal to the saturation value at the seal temperature. Integrating Equation 
(2.6) using the boundary conditions that at r = r b , p = p b and at r = r,-, p = Pi, we get, for a 


parallel seal. 
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For a seal with constant coning slope, the expressions are, 
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In addition to the above equations, the thermodynamic data of the sealed fluid must 
be provided. For the present model, the saturation data, viscosity, density and gas constant 
of the vapor must be provided. These data may be provided as a function subprogram in 
the computer code. The present form of the code uses simple fourth-order spline fits of 
saturation data and incorporates water, nitrogen, oxygen and hydrogen. Additional fluids 
may be easily added and the routine changed without affecting the main code. 


Neglecting convection in the fluid and the inertia terms in the energy equation (2.7), 
the heat conduction rate into the seal plate is a simple expression. 



Since most seals operate at temperatures far below the critical temperature of the 
working fluid, the vapor viscosity is much less than that of the liquid. The heat generation 
rate in the vapor region is therefore much smaller than in the liquid region and can be 
neglected. To evaluate the temperature rise at any radial position r s , we integrate Equation 
(2.28) only over the area of the liquid region of the seal. 
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Assuming angular symmetry, the integral can be expanded to. 
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The seal temperature of interest is the temperature evaluated at the boiling location. 
The pressure at the boiling location must be the same as the saturation pressure of the fluid 
corresponding to the seal temperature there. Recognizing the angular integral as a complete 
elliptical integral of the first kind. Equation (3.1 1) can be simplified to the following, where 
K represents the complete elliptic integral of the first kind, note that in generally the film 
thickness h is a function of the radial position r , 



The torque, Q, required to overcome the viscous friction is usually very small and 
can be found using the following expression, where t 2 q is the shear stress in the 


circumferential direction, 



Using the same argument as before, we can neglect the frictional torque produced in the 
vapor region. We can rewrite Equation (3.13) as, 
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For a parallel seal, the torque required to overcome viscous friction is, 
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For a seal with constant coning slope, the torque required is. 
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(3.16) 


3.4 Method of Solution 

A computer code suitable for PC operation has been developed which analyzes the 
low leakage laminar seals which is documented in Appendix D. The boiling location of a 
low leakage two-phase seal is a operation parameter of significance. When the convection 
terms are negligible, one may calculate the seal temperature based on the boiling location 
alone using Equation (3.1 1). After the seal temperature at the boiling location is found, we 
can determine the liquid and vapor viscosities and densities and the boiling interface 
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pressure, which is the saturation pressure corresponding to the seal temperature at the 
boiling location. 

The solution method is as follows, 

1 . The seal geometry, fluid, inlet and exit conditions, bulk temperature of the sealed 
liquid, seal thermal conductivity and speed are chosen. 

2 . A value of the film thickness (inlet value for the case of a coned seal) is chosen. 

3. Check whether or not the seal is operating at above critical speed. If the seal is, 
calculate where the liquid front is located, determine the seal opening force and then 
go to step 10. If the seal is not, proceed to step 4. 

4. Check whether or not the seal is two-phase. If it is all- liquid, determine the 
pressure profile using Equation (3.1) or (3.3), the seal opening force, and the 
leakage rate using Equation (3.2) or (3.4), and then go to step 10. If the seal is 
two-phase, proceed to step 5. 

5. A boiling interface location is assumed. The seal temperature at the assumed 
boiling location is calculated by numerically integrating Equation (3.12). Note that 
the integrands in Equation (3.12) are singular at the boiling location. An open 
interval integration scheme such as one of the Gauss quadratures performed over 
several sub-intervals, with very small ones near the point of singularity, is 
necessary for an accurate numerical evaluation. 

6. The boiling interface pressure and the fluid properties are calculated from the 
thermodynamic data. 

7. Using the liquid equation, (3.1) or (3.3), and the vapor equation, (3.6) or (3.8), the 
leakage rates are separately determined for both liquid and vapor. These are 
compared and if the vapor leakage rate is higher then the boiling location must be 
shifted towards the inlet and vice versa. 

8 . Steps 5 through 6 are iterated until the proper boiling location is determined. 


76 



9. Pressure profiles are found using Equation (3.2) or (3.4) for the liquid region and 
Equation (3.7) or (3.9) for the vapor region. The opening load is found by 
numerically integrating the pressure over the seal face. 

10. The process is repeated for a different film thickness in order to generate a stiffness 
curve. 
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CHAPTER 4 


TURBULENT TWO-PHASE SEALS - THEORY AND EQUATIONS 
4.1 Introduction 

The full problem of determining the flow through seals in the turbulent regime is a 
difficult one. It involves the inclusion of the nonlinear fluid inertia terms, the consideration 
of fluid friction in which the flow has large velocity components in two directions, and in 
many cases phase change and choking. To make matters worse, the inlet velocities in 
turbulent seal flows are large enough to cause inlet pressure losses to be significant. Since 
no simple boundary condition can be specified at the seal inlet a priori , the problem domain 
must be extended beyond the seal proper, upstream to a location where the fluid properties 
are known. These features of the general turbulent seal problem are in direct contrast to 
those of classical lubrication theory where nonlinear terms are neglected (making analytical 
solutions possible in many cases), where frictional effects are easily handled, and where 
the domain of interest need not extend beyond the flow channel of the device being treated. 
Solution to the full turbulent seal problem requires sophisticated multi-dimensional 
numerical techniques which are difficult to produce and verify and are costly to implement 
on available digital computers. 

Fortunately, for any practical seal, the thickness of the flow channel is very small 
compared to other characteristic dimensions of the device. This fact affords great 
simplification to the model governing equations for the seal proper. Order of magnitude 
arguments show that the component of velocity in the direction through the film is 
essentially zero and that through-the-film variations of temperature and pressure are 
negligible. One may then write the full basic equations of continuity, momentum, and 
energy in terms of shear stresses dropping the terms which are small by the above order of 
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magnitude arguments. It is then possible to integrate through the fluid film to find film- 
averaged properties in a fashion similar to the Karman-Pohlhausen integral treatment for 
boundary layer flows. Shear stresses integrated through the film are resolved into wall 
shear stresses which may be modelled using the film-averaged flow velocities by the semi- 
empirical law-of-the-wall theory. 


By reasonable assumptions yet to be discused, the problem of determining the flow 
from the upstream reservoir to the seal inlet may be similarly simplified. Consideration of 
the extended Bernoulli equation and the First Law of Thermodynamics make it possible to 
relate the velocity, pressure and specific enthalpy at the seal inlet to the known 
thermodynamic fluid state in the reservoir. One can then "jump" from the reservoir to the 
inlet without considering the details of the flow between these stations. This eliminates the 
necessity of extending the problem domain beyond the seal itself and making 
multidimensional numerical calculations there. 

The advantage of this method is obvious, the number of spatial dimensions of the 
problem is reduced by one and the computational domain is restricted to the flow passage 
of the seal. Numerical computation is much simplified and the computer cost to solve an 
individual problem is reduced to a level where seal designers may realistically do parametric 

studies. 

In the past, some investigators have applied integral methods to problems involving 
turbulent flow in narrow channels. Constantinescu and Galetuse [115] and Burton and 
Hsu [116] employed the integral approach in their analyses of turbulent journal bearings. 
In both cases their analyses were correctly done; however, neither group looked beyond the 
journal bearing application for this solution method. In both analyses, it was assumed that 
the gauge pressure on the ends of the bearing were zero. The pressure boundary condition 
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then was independent of the velocity there so that the problem of the development of flow 
from the reservoir to the ends of the bearing was not encountered or discussed by these 
authors. Still others have developed their own methods based on the integral approach and 
used it to analyze a variety of situations. Ng and Pan [1 17] proposed what they called "A 
Linearized Turbulent Lubrication Theory" in which the nonlinear mean fluid inertia terms 
were dropped from the film-averaged turbulent equations of motion. Linearized turbulence 
is in this sense logically inconsistent; if the flow is turbulent, the Reynolds number is large 
enough that the mean fluid inertia terms play a significant role. Hirs [118] developed yet 
another approach to turbulent lubrication flow. His bulk flow equations are derived from 
the film-averaged turbulent equations of motion. The turbulent wall shear stresses are 
related to the mean velocities by a correlation made by considering experiments with 
turbulent flow in narrow channels. Although the mean fluid inertia terms are present in the 
equations, Hirs’ formulation mixes these terms with those associated with the wall shear 
stress terms making the role of the inertia terms in any given problem difficult to assess. 
Extension of Hirs’ method to situations which require the inclusion of the energy equation 
would be extremely awkward. 

The underlying philosophy of treating turbulent lubrication in this work is to rely on 
the fundamental equations. Wall shear stresses are to be calculated by friction factors 
which have been used extensively and shown to enjoy a wide range of applicability. All 
the important terms are retained and the inertia terms are isolated from the other terms in the 
equations. Many problems in turbulent lubrication theory may then be approached in a 
simple and consistent way by this method. 
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4.2 Derivation of the Basic Equations for Turbulent Seal Flow 


The integral approach will now be applied to the basic governing equations for both 
the radial face seal and annular seal geometries. In these geometries, a homogeneous fluid 
is assumed to flow. In the case of two-phase flow, the mixture of liquid and vapor is taken 
as a pseudo-fluid with properties which are averages of the properties of the component 
phases. Turbulent flow through these seals is assumed to be adiabatic. This assumption 
does not mean that no heat transfer takes place between the fluid film and the seal itself, but 
rather that the amount of heat transferred from the fluid to the seal over some specified time 
interval is small compared to the amount of heat convected by the flowing fluid or 
generated in the film during that same time interval. The theoretical adiabatic limit is closely 
approached in seals over the entire turbulent range where the leakage rates are large. 
Products of average quantities will be taken as equal to the average of their products in the 
convective acceleration terms. This assumption is not restrictive in turbulent flows where 
the velocity profiles for pressure-driven flow are very blunt 

Basic Governing Equations for Radial Face Seals 

A schematic representation of a radial face seal is shown in Figure 4. The seal is 
composed of two rings, one stationary ring mounted to the fluid barrier known as the seal 
seat, and one ring mounted to the shaft and rotating with it known as the primary ring. 
Fluid in the clearance space between the rings may velocity components in both the 
circumferential and radial directions. Radial flow is pressure driven, while circumferential 
flow is driven by shear from the rotation of the shaft. Depending on the application, one of 
two different seal configurations may be used. If high pressure is felt on the inside radii of 
the seal rings, the configuration is referred to as an inside seal and the radial flow is 
outward; if high pressure is felt on the outside radii of the seal rings, the seal configuration 
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is referred to as an outside seal and the radial flow is inward. The outside seal 
configuration is the more commonly used one and that which is depicted in Figure 4. 




Figure 4-1: Face seal geometry 

(not to scale). 


The opposing surfaces on the face seal rings on either side of the clearance space need 
not be parallel. It will be shown later that it is advantageous to taper the moving ring seal 
face as seen in the schematic diagram. This tapering is generally referred to as coning by 

seal designers. 


By restricting the relative motion of the seal rings to be along their common axis only, 
the problem becomes much simplified; the geometry is axisymmetric and terms for angular 
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variations in the equations ate dropped. Further, time derivatives in the basic equations 
will be dropped as weU. The continuity equation may then be written as: 

7 57< r P u) + S <pv * ) = 0 

where v z is the axial velocity in the z direction. 

Upon integrating from the seal seat to the primary ring, the contribution from the 
second term is zero since the walls are impermeable, making the axial velocity vanish at 
both surfaces. Because the film height changes with the radial position, care must be taken 
to evaluate the remaining integral. Applying Leibnitz’s rule, this integral may be expressed 

as: 

r h(r) . i a r h(r) , ldh 

I 75 f( r pu)dz = ?5F I (r pu) dz - r pul h(r) r ^ 

J o 

where the last term is zero because of the no-slip condition. In the face seal geometry, the 
radius is never zero so that this equation may be integrated to give: 

m = 27trhG (4-1) 

■he film thickness, h, being a known function of radial position. In the above equation, the 
radial velocity, and hence the mass leakage rate and mass velocity, may take either a 
positive or negative sign. In an inside seal the fluid flows radially outward giving a 
positive sign, while in an outside seal the fluid flows radially inward giving a negative 
sign. This convention is used throughout the face seal analysis. 

The radial momentum equation can be processed in a similar way. The general form, 
dropping the angular variations and terms containing the axial velocity, is. 
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(4-2) 
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The axial velocity is necessarily zero at the walls since they are impermeable. It is 
presumed that the axial velocity between the walls can not become so large as to be 
comparable with the velocity components in the other directions. The stress component o n 
is dominated by the pressure. This stress component will be taken to be equal to the 
negative of the pressure. The stress component will be small compared to the other 
components since the mean velocity gradients in the circumferential direction are small 
compared to those in the other directions. 


Special care must be taken to integrate the centrifugal acceleration term through the 
film. In turbulent Couette flow, the velocity profile for the circumferential velocity is "S" 
shaped. If the "S" is shallow, the circumferential velocity can be approximated by w(r,z) = 
torz/h. This expression is then substituted into the centrifugal acceleration term of equation 
(4-2) and the integration carried out Using the stated assumptions and integrating gives: 
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= 0 


(4-3) 


where is the wall shear in the radial direction. The positive sign on the shear term 
applies to outward radial flow (inside seal), while the negative sign applies to inward radial 
flow (outside seal). 


The energy equation can be written for a ring-shaped control volume of width Ar 
which occupies the entire space between the seal set and primary ring as: 

W s = - m(e l r + A r-el r ) 
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where e is the specific energy of the fluid, or: 


2 2 

e = i + u /2 + w /2. 

The left hand member of the above equation represents work down to the fluid by shear or: 
W s = 2 it ( A r) cort z0 


Substituting, and taking the limit as the ring width approaches zero gives the adiabatic 


energy equation for the face seal: 
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(4-4) 


where lh is the wall shear stress in the circumferential direction evaluated at the shaft 
surface. The equations (4-1), (4-3), and (4-4) comprise the set of basic descriptive 
equations for the axisymmetric radial face seal. The circumferential equation of motion for 
this situation is trivially satisfied. 


Basic Equations for Turbulent Flow in Annular Seals 

A schematic representation of an annular seal is shown in Figure 4-4. This seal is 
made up of the shaft itself and a close fitting opening around it. Fluid in the clearance 
space is forced through the seal axially by pressure and is dragged along in the 
circumferential direction by the shearing action of the shaft 
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Figure 4-2: Annular seal geometry (not to scale). 


The through-the-film integration for the annular seal geometry is proceeded in a 
similar fashion to that for the radial face seal and all the same assumptions and mathematical 
techniques applied to the previous case will again be applied. However, the film-averaged 
annular seal equations given here will be more general and will include circumferential 
variations. Because the film dimension in the radial direction is very small, the centrifugal 
acceleration can not build up any appreciable pressure difference. Centrifugal acceleration 
effects in this case will be neglected. 

Since the film thickness is very small compared to the shaft radius, curvature effects 
can be neglected and the problem domain may be unwrapped; the coordinate system 
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reduces to a Cartesian system. Integrating through the film of the continuity equation and 
simplifying gives (Note that here we use u and w for axial and circumferential velocity, 

respectively): 

^-(puh) + ^(pwh) = 0 (4-5) 

while the equations of motion in the axial and circumferential directions are, respectively. 

puh^ + P Wh ^ 1^" = ' ( X rzlR+h “ T izIr) ~ (4-6) 


and 
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Finally, the adiabatic energy equation becomes: 


x r 0 | R 0) R 

= puh (i + u 2 / 2 + w / 2 ] + pwh (i + u / 2 + w 72 ) 


(4-8) 


The equations (4-5), (4-6), (4-7), and (4-8) make up the set of descriptive equations 
for the flow within the annular seal. In all the above equations for an annular seal, the 
variation in radius is neglected. The mean radial position within the film is approximated 
by the shaft radius. 
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4.3 Determination of Wall Shear Stresses 


Many different approaches to finding wall shear stresses in turbulent flow have been 
proposed. Burton [119] gives a good catalogue of these approaches, summarizing their 
origins and salient features. Of all these methods, the semi-empirical law-of-the-wall 
theory is used here to give wall shear stresses in terms of the components of the mean 
channel velocity. As stated previously, this choice of turbulent friction theory has 
advantages over other theories, namely, the theory enjoys wide applicability and the mean 
inertia terms are separated from other terms in the equations of motion. 

From variations in fluid density or channel height, the mean flow velocities and hence 
the wall shear stresses may change over the area of the seal. These variations may be 
significant and must be dealt with. In all practical cases, the flow channel height is small 
compared to other dimensions of the seal and that that channel height varies slowly along 
the flow path. The flow may then be taken to be quasi-fully developed, and the law-of-the- 
wall theory, applicable to fully-developed flow, may then be applied locally to give the wall 
shear stresses in terms of the fluid properties there. The relation: 

f ,,2 
X= 8 pV 

then holds, where t is the local wall shear stress, V is the mean velocity in the direction in 
which the shear acts, and f is a friction factor. Correlations exist which specify the friction 
factor, depending on the flow situation encountered. 

Under typical operating conditions, the flow could have large mean velocity 
components in two directions simultaneously, giving a complex three-dimensional 
turbulent structure in which coupling does exist between the velocities in these different 
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directions. Various treatments for flows of this type have been used. For example. Black 
and Jenssen [120] apply an elaborate correlation proposed by Yamada [121] to analyze 
turbulent annular seals. White [122] lists several approaches which account for turbulent 
flow coupling. However, he also comments that, at present, no good theory exists for 
treating flows of the type just described. In light of this fact, the use of simple turbulence 
models which give reasonable results are warranted. 

One particularly simple model proposed for coupled turbulent flow treats the flows 
associated with the different velocity components separately [122]. This approach greatly 
simplifies the calculations. Trial studies indicate that the fluid property profiles and leakage 
rates are rather insensitive to variations in the friction factor values m either direction, 
suggesting that the uncoupled flow model is indeed adequate. 

The channel geometry for seals is always such that the aspect ratio, or ratio of channel 
height to width, is extremely large so that the simple hydraulic radius concept is not 
applicable. Instead of using a hydraulic radius, it is more accurate to model the channel as 
one that is infinitely wide. Special friction factor correlations have been made for 
hydraulically smooth-two-dimensional channels. For pressure-driven or Poiseuille flow, 
the proper friction factor correlation in turbulent flow is reported by White [122] to be. 

-L = 0.8839 In (Re h Vf) + 0.142 ( 4 ‘ 9 ) 

Vf 

The Reynolds' number based on film height is defined as: Re h = . For shear-driven, or 

Couette, flows the mean velocity profile across the channel assumes an "S" shape. The 
law-of-the-wall under a zero prevailing pressure gradient may be manipulated to give the 
following friction factor correlation as derived by that same author [122]: 
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-L = 1.768 In (Re h Vf) + 0.831 

Vf 


(4-10) 


where action, i, has been shown that this correlation is still valid if the 

prevailing pressure gradient is negative or tnoderatively positive. Solution of these 
transcendental equations for the friction factor may be easily accomplished numerically by 
Newton's root-finding method since the local velocity and fluid properties are known. In 
the two-phase regime, the mixture density is used and the viscosity is taken to be the 
volume-weighted average of the viscosities of the component phases. 

4.4 Calculation of Inlet Thermodynamic State 

It is assumed that the upstream reservoir area is large compared to the flow area at the 
seal inlet so that the velocity of the reservoir fluid is essentially aero to an observer in the 
laboratory. Temperature and pressure in this body of fluid far from the seal are measurable 
and uniform. The possibility does exist that parts of the seal at an elevated temperature can 
transfer heat to the reservoir fluid, but this effect is assumed to be negligible. The zone of 
influence of hot seal elements on the temperature of the incoming fluid will not extend far 
upstream into the reservoir. The fluid entering the zone where heat transfer from the 
mechanical paris of die seal could be appreciable is travelling a. a velocity which is a large 
fraction of the inlet velocity. It is then reasonable to assume that the flow is nearly 
adiabatic, as was done for flow in the body of the seal. 

With the above assumptions, the first law of thermodynamics may be used to give: 

2 
« i 

i i= I oo- (4-11) 

since the difference in potential energy between these two stations is very small. Again 
neglecting any small potential energy change, the extended Bernoulli equation may be 
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written for a streamline from a point far upstream in the reservoir to just inside the seal 
inlet: 

(4-12) 

oo ^ 

where k is an inlet loss coefficient which accounts for the irreversible pressure loss due to 
sudden contraction of the flow area. The loss coefficient multiplies the mean kinetic energy 
per unit mass at the inlet. Values for this coefficient may range from 0.01 for a well 
rounded inlet to 0.5 for a square-edged inlet. 

Nowhere in the above equations has the circumferential, or swirl, velocity appeared. 
It has been assumed that the work done by rotating seal elements on the incoming fluid by 
shear is approximately equal to the kinetic energy increase of that fluid by increased swirl 
velocity. These terms then cancel each other and so do not appear in the governing 
equations. This statement is equivalent to assuming that there is a negligible pressure 
difference through thin boundary layers surrounding rotating seal parts. 
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Figure 4-3: Typical trajectory of a fluid element 
S travelling through the seal in the T-s plane. 


If the reservoir fluid is very close to its saturation state, it is certainly possible for 
vaporization to occur as the fluid travels from the reservoir to the inlet as the pressure tt 
experiences decreases. Such a case is shown on the T-s diagram for a typical substance in 
Figure 4-3 where the first segment of the curve indicates the change of state from the 
reservoir to the inlet and the second segment shows the change in state through the seal 
proper. The integral that is the left hand member of the extended Bernoulli equation ts then 
difficult to evaluate because it is path-dependent. Fortunately, for most design situations of 
interest, the amount of vaporization from inlet -pressure drop will be small; most of the 
boiling of liquid will occur within the body of the seal. Taking the fluid quality just inside 
the inlet to be small, the integral of Pdv, or the work done on a system of fluid travelling 
from the reservoir to the seal inlet may be approximated by: 
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in other words, the path taken through p-v space is taken as a straight line. Integrating 
equation (4-12) by parts, and substituting from equation (4-13) gives the following: 
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(4-14) 


This equation, along with the First Law, from a set of equations which specify the 
thermodynamic state of the fluid at the inlet 

Simultaneous solution of equations (4-11) and (4-14) is accomplished by iteration. 
The calculation is begun assuming that the fluid remains a liquid through the inlet process. 
The specific volume then remains a constant over a streamline and the inlet velocity is 
found from continuity for a presumed mass flow rate. (To avoid confusion, it should be 
noted that the mass flow rate is determined independently from this inlet state calculation 
from pressure boundary conditions imposed upon the entire flow problem. Solution 
methods to find this mass flow rate will be outlined in subsequent chapters.) Since the 
specific volume is known for the entire inlet pressure drop, equations (4-11) and (4-14) 
may be used separately to determine the specific enthalpy and pressure at the inlet. These 
two fluid properties fully determine its thermodynamic state. If this calculated state indeed 
corresponds to a subcooled liquid, the initial premise was correct and the inlet state is now 
known. If the calculated inlet state corresponds to a saturated mixture, however, the 
iteration must continue. To find the inlet state when that state is under the saturation dome, 
the estimate for the inlet state is used to find the inlet specific volume. From this estimate 
and the presumed mass flow rate, an updated inlet velocity is calculated. Substituting the 
updated inlet velocity and specific volume into equations (4-11) and (4-14) gives updated 
values of pressure and enthalpy which define a new updated thermodynamic state. This 
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new estimated state is compared to the previously calculated state. Calculation continues 
until the relative difference between successive inlet pressure and enthalpy estimates falls 
below some tolerable error bound. By the nature of the relation between the specific 
volume and velocity, this procedure will always converge to a solution. 
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CHAPTER 5 


LEAKAGE AND AXIAL STIFFNESS 
IN TURBULENT RADIAL FACE SEALS 

5.1 Introduction 

In this chapter, turbulent radial face seals will be analyzed using the basic equations 
derived in section 4-2 to describe the flow in the seal proper and other equations derived in 
section 4-5 to relate the inlet conditions to the reservoir conditions. An idealized seal will 
be considered in which the seal rings: 1) have hydraulically smooth surfaces, 2) are 
constrained to move only in the axial direction relative to each other rendering the geometry 
axisymmetric, and 3) move so slowly that squeeze-film damping effects are negligible. In 
reality, the seal rings in most applications are polished by methods like those used to 
manufacture telescope mirrors and do approach the hydraulically smooth limit for 
reasonable clearance gaps. The seal rings are connected to nearby machine elements by 
elastomeric materials and so actually execute relative motions in directions perpendicular to 
the axis of the seal; however, these motions are dominated by the axial motion. Finally, it 
should be noted that damping, though quite important in the dynamics of face seals, does 
not determine the conditions for seal collapse. One need not be concerned with the degree 
of damping to investigate the axial stability of the seal. Thus, the highly idealized model 
presented is not restrictive and gives much useful information about the performance of real 

face seals. 

One may generate a curve of opening force as a function of mean film thickness for 
any set of operating conditions. Assuming that the seal has sought a film thickness where 
the forces on the primary ring are in equilibrium, the slope of the opening force curve at 
that point, multiplied by minus one, gives the restoring force due to a small displacement 
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from the equilibrium position. The entire response of pressure in the fluid film can then be 
modelled to the action of an equivalent single axial spring with a variable, but known, 
spring constant. It is possible for that spring constant to become negative under some 
operating conditions. If the equivalent spring constant does become negative, the net 
pressure on the primary ring acts not to restore the seal to equilibrium but instead acts to 
collapse the seal. 

It will be shown that seal leakage and opening force vary with the rotation speed of 
the spinning face and with the thermodynamic state of the incoming fluid in the two-phase 
regime. In many situations, the inlet drop is a significant fraction of the total pressure drop 
through the seal. Further, it will be seen that the inlet drop can give rise to unexpected seal 
behavior. 

Choking may occur in the seal if the back pressure is below a critical value. At 
choking conditions, the flow velocity is equal to a sonic velocity just at the seal exit and the 
exit pressure no longer matches the back pressure. Pressure trajectories through the seal 
are shown in Figure 5- 1 for choked and unchoked flow. It is possible for the fluid at the 
seal exit to be vapor, two-phase mixture, or even pure liquid under choked conditions. The 
all-liquid choked case may occur when the rate of viscous dissipation is small and the 
degree of inlet subcooling is large. This mode of choked flow behaves much differently 
than the two-phase or all vapor choked flow. The previous analysis of Hughes and Beeler 
[89] allows only for saturated liquid entering the seal and so can not predict this unusual 
mode of choking. The later paper by Beatty and Hughes connects these deficiencies [90]. 
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Figure 5-1: Typical pressure profiles through a face seal. 

Curves "a" and "b" show unchoked flow with exit and back 
pressures equal. Curve "c" shows the choked flow line 
and expansion wave with the back pressure lower than the 
critical pressure. Curve "d" shows the cavitation line 
for an inside seal. 


As discussed by Osterle and Hughes [103], it is possible for cavitation to occur in an 
inside seal if the quality is low throughout the flow passage and the rotation speed is 
sufficiently high. In this case the pressure drops to the value of the back pressure within 
the interior of the seal, rupturing the film without boiling. The net opening force under 
such circumstances would be extremely low. This possibility is depicted in Figure 5-2. 

The criterion for cavitation is that the pressure gradient goes to zero as the pressure 
approaches the back pressure at some point within the interior of the seal. Although this 
condition is not complicated, determining whether cavitation occurs is very difficult in that 
it depends on the specific parameters of the problem. If sufficient boiling occurs within the 
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seal, the fluid density is slow so that the centrifugal inertia term is small and cavitation does 
not occur. No attempt will be made to do a systematic study of cavitation in turbulent seals 
in this analysis. 

5-2 Specialized Equations for Liquid, Two-Phase, and Vapor Flow 

The basic equations are now adapted for liquid, two-phase, or vapor flow. As stated 
previously, two-phase mixtures are treated as homogeneous fluids. The continuity 
equation does not change with change of phase and is: 

m = 27trhG (5-1) 

Remember that the mass leakage velocity can be positive or negative. It is taken as positive 
in an inside seal and negative for an outside seal. 

Liquid Flow 


If the fluid is entirely liquid, it can be assumed with little loss of accuracy that the 
density is constant with a value corresponding to the density of saturated liquid at the same 
temperature. This affords some simplification of the radial momentum equation which then 
can be written: 
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where the shear stress term takes a negative sign for flow radially outward. Further, the 
energy equation becomes: 
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Note tot the derivative of the Him thickness with radius in these equations is simply the 
tangent of the coning angle, (3. 


Since the liquid compressibility is neglected, the sonic speed for the liquid is infimte 
and there is no singularity in either to pressure gntdient or to enthalpy gradient. Normal 
choking, where to fluid velocity reaches to acoustic wave speed in to medium, then can 
not be predicted from to equations above. I. should be noted to. to acoustic wave 
speeds for most liquids are so high tot choking in a mode as described above is very 

seldom experienced. 


Liquid - Vapor Flow 


For a homogeneous saturated two-phase mixture in thermodynamic equilibrium, the 
specific enthalpy and specific volume may be written in teims of quality: 

— = V = Vf + Xv fg 

p 


and. 


i = i f + Xi fg 


Further, the radial derivatives of these may be rewritten in tetms of other quantities by 
application of the chain rule, for example: 


di_ 

dr 


di 

dp 


U 1 t-rr I 

- = dr 


dp . dX, 
- + ^ dF 


where. 
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di 

¥ 


'k = 




The radial velocity in the radial equation of motion and the energy equation may be 
expressed in terms of the mass leakage velocity, G. Using this and the relations among the 
saturation properties above, these equations may be solved for the radial derivative of 
pressure explicitly, assuming G is known: 


dp = 1 OH 

dr i . § \ r \ 


L dh 

+ h dr 




+ G 


2Vfg 

*fg 


tor 

Gh Xze + 


2 

2wr ± 
3 



(5-4) 


where 


<j> = - G 


2 dv . 



di 

'3p 


Again, the radial derivative of film thickness is simply the tangent of the coning angle (3. 
In the above expression the negative signs are taken for the radial shear stress terms if the 
fluid is flowing radially outward. Once the pressure gradient is known, the mixture quality 
gradient may be expressed as: 


2 

dX _ 1 tor _ 2to r 

dF " ' ifgl" Gh ze 3 


di 

'3p 


2v 


dp 



(5-5) 
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where, once again, the negative sign applies to the radial shear stress term for radial 
outward flow. Since temperature is only a function of pressure at saturation, the following 
is true: 


di_ - fe + jl + i — 

dr (dp dp / dr fg dr 


(5-6) 


where the pressure gradient is known from (5-4) and the quality gradient is known from 
(5-5). Note that equations (5-4) and (5-6) give finite gradients, as in the all-liquid case, 
unless <(> becomes greater than or equal to unity. The quantity <t> corresponds to the square 
of the Mach number in an ordinary gas dynamics problem; later, it will be related to a 
characteristic wave speed in a two-phase medium. 

All-Vapor Flow 


If one assumes an ideal gas with constant specific heats, the continuity, axial 
momentum, energy, and state equations may be combined in a simple way to give the 
gradients of the pressure and enthalpy. The pressure gradient may be found explicitly as: 

dp _ l_/lpM I/ r dhl 

* i-M 2 \ r ' h *J 


4 - 



(l + 2 (y- 1)M 2 ) 


(5-7) 


±i(l + (T-DM 2 )x„- <T ~o h V ^— 


and, the enthalpy gradient, in terms of the pressure gradient expression above, is: 
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(5-8) 
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— 1 T X J + 



with the radial derivative of film thickness specified as before. Here, M is the Mach 
number. 


5-3 Solution Procedure 


Once a mass flow rate is specified, the appropriate set of governing equations may be 
integrated to give pressure and enthalpy at all points interior to the seal since the upstream 
reservoir conditions are known and definite relationships exist between the conditions of 
the reservoir and of the seal inlet. Remember that the mass velocity, G, is related to the 
flow rate and the local conditions through equation (5-1). These equations are nonlinear 
and so can not be treated analytically; they are instead integrated in this study numerically 
by a standard fourth-order Runge-Kutta technique. The choice of the set of governing 
equations used is made by examining the thermodynamic state of the fluid at the last node 
where the solution is known. Equations (5-2) and (5-3) apply to the subcooled liquid, 
equations (5-4) and (5-6) apply to a saturated mixture, and equation (5-7) and (5-8) apply 
to a superheated vapor. 


A solution to the problem as posed is found when a specified mass flow rate either 
gives an exit pressure equal to the back pressure, or is at its maximum (choked flow). 
Referring back to Figure 5-1, these possibilities are shown as curves "a", "b", and "c". 
Since the critical back pressure and mass flow rate are not known a priori , a test for 
choked flow at the outset of the solution procedure is necessary. 

When the flow is choked, the radial derivative of pressure is infinite at the seal exit. 
This condition corresponds to having the denominator of equation (5-4) for two-phase flow 
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or the denominator of equation (5-7) for vapor flow go to zero just at the exit. The choking 
condition may be found numerically within preset error limits by a method similar to the 
bisection method for finding roots to algebraic equations. First, a trial mass flow rate is 
chosen to insure that the denominator of the appropriate pressure gradient expression does 
not go to zero or change sign at any node point in the seal. Next, another larger trial mass 
flow rate is chosen to insure this sign change in the denominator takes place at some 
interior node. The critical mass flow rate then must occur within the interval defined by the 
mass flow rates discussed above. The original interval may then be subdivided and the 
denominator of the pressure gradient expression tested at all nodes down the axis for the 
intermediate mass velocity. This will further narrow the critical flow rate to one of the 
subintervals. This process is continued until the possible error in the critical flow rate is 
within tolerable bounds. The critical back pressure is the exit pressure for the critical mass 

flow rate. 

If the given back pressure is flower than the critical back pressure found by the 
method described, the problem is solved; the mass flow rate is at its critical value and the 
procedure is terminated. If not, the problem takes the form of the standard shooting 
method for a two-point boundary value problem, with mass low rate as the shooting 
parameter. Finally, the leakage rate is reported and the opening force is calculated as the 

are integral of the seal pressure. 

5-4 Discussion of All-Liquid Choking 

It is possible to have an anomalous situation where the flow is choked but remains 
liquid up to the seal exit. This is likely to occur if the rate of viscous heat generation is 
small and the degree of subcooling at the upstream reservoir is large. The liquid pressure 
will fall, at the exit, to the saturation pressure for the local temperature, but no lower. 
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Beyond the exit, the liquid flashes immediately into vapor, an outside observer would not 
see any liquid issuing from the seal. 

This peculiar choking mechanism may be understood through the careful study of 
longitudinal wave propagation in saturated two-phase mixtures as discussed by Pai [126], 
In this instance, there are two modes of wave propagation in the medium. One mode is that 
of ordinary sound waves travelling through the vapor phase, while the other, the so-called 
vaporization wave mode, is peculiar to the saturated mixture. Unlike the ordinary sound 
wave, the vaporization wave brings about a change of phase. One then would expect that 
the propagation speed of the vaporization wave is entirely different than that of an ordinary 
sound wave. In fact, the speed of the vaporization wave is given by the square root of the 
expression <> in (5-4). Note that this quantity has a finite limit, even as the mixture quality 
approaches zero. Generally, the vaporization wave speed is much lower than that of the 
ordinary sound wave under the same background conditions. Thus, the speed at which 
information can be relayed in a saturated mixture is limited by the vaporization wave speed. 

Anomalous all-liquid choking occurs when very slightly subcooled liquid is moving 
faster than the local vaporization wave speed at the exit. If one imagined that the flow rate 
could increase, the exit pressure would have to decrease in response and the liquid at the 
exit would become saturated. But, the saturated mixture exit velocity is limited by the 
vaporization wave speed. Since the vaporization wave speed is smaller than the presumed 
exit velocity for the mixture, this situation is physically impossible. Hence the flow in this 
instance is choked, but is subcooled liquid throughout the seal. 
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5-5 Numerical Examples and Discussion 


Numerical examples are presented based on two applications, namely: a nominal 
design for the interstage seal of the Space Shuttle Main Engine High Pressure Oxidizer 
Turbopump, and a high pressure water pump. In both cases, the seal inner radius is taken 
to be 0.0430 meters, while the outer radius is taken to be 0.0469 meters. The flow inlet is 
taken to be square-edged with a nominal inlet loss coefficient of 0.5 in all cases. 

Figures 5-2 through 5-5 show fluid property profiles for the nominal Space Shuttle 
interstage seal under various possible conditions. These given conditions do not 
necessarily correspond to typical operating values, but rather demonstrate some extremes of 
face seal behavior. In these examples, cryogenic oxygen is taken as the sealed fluid. 

Figure 5-2 shows an unchoked flow in an outside seal configuration. The flow 
enters as a two-phase mixture and remains so throughout. The pressure profile is smooth 
with a finite gradient everywhere. The predicted leakage in this case is 3.28 x 10*2 kg/s 
and the predicted opening force is 3670 N. 

Figure 5-3 shows profiles for an outside seal with reservoir conditions identical to the 
previous case, but a back pressure below the critical back pressure. The pressure gradient 
is extremely steep at the seal exit, demonstrating two-phase choked flow. The flow rate is 
maximum for the given reservoir conditions with a value of 3.36 x 10' 2 kg/s. The opening 
force is 3600 N, less than the previous case. This is expected since the pressure at any 
point is less than the corresponding point in the unchoked case. 

Figure 5-4 shows nearly straight line profiles for the all-liquid choked flow case for 
an outside seal configuration. Although the fluid exit velocity is much lower than the sonic 


106 



speed of the liquid, the mass flow is at its maximum value and is independent of back 
pressure, provided it is below a critical value as discussed previously. The leakage rate 
here is 7.25 x 10 -2 kg/s, nearly twice that of the previous cases. The opening force is 3040 
N. Large inlet losses partially account for this low opening force value. 

Figure 5-5 shows an unchoked flow in an inside seal configuration without 
cavitation. The operating conditions are identical with those of the case shown in Figure 5- 
2. One notes that the leakage rate is 3.85 x 10 2 kg/s, significantly higher than in the 
outside seal configuration, while the opening force is 3620 N, a value less than that of the 
outside seal case. These differences are a manifestation of the centrifugal inertia. The 
effect of the centrifugal inertia is to drive, rather than impede, the flow to the exit giving a 
higher leakage rate. The higher leakage, in turn, gives larger inlet losses so that the 
pressure throughout the seal is lower than in the outside configuration. 

The examples shown in Figures 5-6 through 5-16 are for cryogenic oxygen flow in 
an outside seal. In each case, the seal leaks fluid to a vacuum for choked flow. The final 
examples shown in Figures 5-17 and 5-18 are for high pressure water leaking into a space 
where the ambient pressure is atmospheric. 
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Figure 5-2: Fluid property profiles for an unchoked outside seal. 

Sealed fluid is Oxygen 
Film Thickness, h = 1.00 x 10‘^ m 
Shaft Speed, to = 30,000 RPM 
Reservoir Temperature, Too = 150.0 K 
Reservoir Pressure, Poo = 4.30 MPa 
Back Pressure, Pb = 2.00 MPa 

Coning Angle, P = 0.0 radians 
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Figure 5-3: Fluid property profiles for a choked outside seal. 

Sealed fluid is Oxygen 
Film Thickness, h = 1.00 x 10' 5 m 
Shaft Speed, co = 30,000 RPM 
Reservoir Temperature, T« = 150.0 K 
Reservoir Pressure, Poo = 4.30 MPa 
Back Pressure, Pb = 0.0 MPa 

Coning Angle, P = 0.0 radians 
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Figure 5-4: Fluid property profiles for an outside seal. 

with all-liquid choked flow. 

Sealed fluid is Oxygen 
Film Thickness, h = 1.00 x 10' 5 m 
Shaft Speed, to = 15,000 RPM 
Reservoir Temperature, Too = 125.0 K 
Reservoir Pressure, Poo = 4.30 MPa 
Back Pressure, Pb = 0.0 MPa 
Coning Angle, P = 0.0 radians 
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Figure 5-5: Fluid property profiles for an unchoked inside 

seal without cavitation. 

Sealed fluid is Oxygen 
Film Thickness, h = 1.00 x 10"^ m 
Shaft Speed, 0) = 30,000 RPM 
Reservoir Temperature, Too = 150.0 K 
Reservoir Pressure, Poo = 4.30 MPa 
Back Pressure, Pb = 2.00 MPa 

Coning Angle, P = 0.0 radians 
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Figure 5-6 shows the effects of rotation speed and film thickness on the leakage rate. 
As the speed is increased or the mean Film thickness is reduced, more vapor is generated in 
the seal. Choking is promoted and the leakage is reduced. In addition, for an outside seal 
configuration, increasing the rotational speed increases the contribution of centrifugal 
inertia which retards the flow and so further reduces leakage. Similar results are obtained 
for all cases where the degree of subcooling in the upstream reservoir is moderate. 

Figures 5-7 through 5-11 show the effects of film thickness and rotational speed on 
opening force for various values of initial subcooling. The interesting features of these 
curves can be largely explained by the competition of two opposing effects, the decrease of 
opening force by dissipation and the increase of the inlet pressure level through reduced 
inlet losses. In rough terms, dissipation and pressure loss promote boiling giving 
increased quality. If the mixture quality is increased, its density decreases. For a nearly 
constant area channel, the convective acceleration of the fluid is likely to be larger and the 
wall shear stresses increase. The pressure gradient must then be larger to offset these 
increases. Larger pressure gradients drive the pressure profile down so that the area under 
the pressure-radial position curve is reduced. The opening force is directly proportional to 
the area under that curve. On the other hand, boiling reduces the leakage rate and so 
decreases the inlet flow velocity. The inlet losses are reduced so that the inlet pressure level 
increases. This increase raises the pressure profile and so increases the opening force. 
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Figure 5-6: Effects of rotation rate and film thickness on 
leakage with 0.3 K subcooling. 
Sealed fluid is Oxygen 
Reservoir Temperature, To© — 150.0 K 
Reservoir Pressure, P<» = 4.3 MPa 
Back Pressure, Pb = 0.0 MPa 

Coning Angle, P = 0.0 radians 
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Fieure 5-7: Effects of rotation rate and film thickness on 
leakage with 0.3 K subcooling. 
Sealed fluid is Oxygen 
Reservoir Temperature, Too = 150.0 K 
Reservoir Pressure, Poo = 4.3 MPa 
Back Pressure, Pb = 0.0 MPa 

Coning Angle, P = 0.0 radians 
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Figure 5-8: Effects of rotation rate and film thickness on 
” opening force with 3.3 K subcooling. 

Sealed fluid is Oxygen 
Reservoir Temperature, Too = 147.0 K 
Reservoir Pressure, Poo = 4.3 MPa 
Back Pressure, Pb = 0.0 MPa 
Coning Angle, (5 = 0.0 radians 
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Figure 5-7 gives stiffness (i.e., force-displacement) curves for a slight initial 
subcooling. At large film thicknesses, the three lowest speeds give opening forces that 
increase with speed. This happens because the inlet loss reduction outweighs the force - 
reducing effects of dissipation. At small film thicknesses, inlet loss effects are small since 
the leakage rates are small. Then, the opposing, force-reducing effect dominates and the 
relative order of the opening forces with speed are reversed. Under the highest speed 
condition the leakage rates are yet smaller than the cases previous, so the force-reducing 
effect dominates the inlet loss effect even at large mean film thicknesses. The high speed 
stiffness curve then lies below the others over most values of film thickness. Note that the 
curves give positive stiffness (i.e., a negative slope for the force-displacement curve at the 
point of interest) for large film thicknesses, but give neutral or negative stiffness for small 

thicknesses. 

Figure 5-8 gives stiffness curves for an initial subcooling somewhat greater than the 
case before it. Results are similar to those of the previous case, but the inlet loss effects are 
relatively less strong. The curves exhibit stiffnesses that are everywhere less positive than 

before. 

Figure 5-11 gives stiffness curves for a still greater initial subcooling. The inlet loss 
effect is almost completely dominated. The stiffness curves are flat or slope upward, 
indicating seal collapse is imminent. 

Figure 5-10 gives stiffness curves for yet a greater initial subcooling. The increase in 
leakage rate due to the suppression of the boiling is large enough to cause the inlet loss 
effect to again be dominant. This plot is similar to that shown in Figure 5-7. 
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Figure 5-9: Effects of rotation rate and film thickness on 
opening force with 5.3 K subcooling. 

Sealed fluid is Oxygen. 
Reservoir Temperature, Too = 145.0 K 
Reservoir Pressure, Poo = 4.3 MPa 
Back Pressure, Pb = 0.0 MPa 

Coning Angle, |3 = 0.0 radians 
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Figure 5-10: Effects of rotation rate and film thickness on 
opening force with 10.3 K subcooling. 

Sealed fluid is Oxygen. 
Reservoir Temperature, T„ = 140.0 K 
Reservoir Pressure, Poo =4.3 MPa 
Back Pressure, P 5 = 0.0 MPa 

Coning Angle, P = 0.0 radians 
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Figure 5-11: Effects of rotation rate and film thickness on 
leakage rate with 50.3 K subcooling. 

Sealed fluid is Oxygen. 
Reservoir Temperature, Too = 100.0 K 
Reservoir Pressure, Poo = 4.3 MPa 
Back Pressure, Pb = 0.0 MPa 

Coning Angle, P = 0.0 radians 
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Figure 5-12: Effects of rotation rate and film thickness on 
opening force with 50.3 K subcooling. 

Sealed fluid is Oxygen. 
Reservoir Temperature, Too = 100.0 K 
Reservoir Pressure, Poo = 4.3 MPa 
Back Pressure, Pb = 0.0 MPa 

Coning Angle, (3 = 0.0 radians 
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Figures 5-11 and 5-12 give leakage rate and stiffness information for a very high 
degree of initial subcooling. Whereas in the previous cases the flow was in the two-phase 
regime upon entering the seal or entered as a liquid but became two-phase soon after, in 
these plots the flow remains liquid throughout much of the flow passage. In fact, the flow 
is all liquid for all film thicknesses shown for the lowest three speed conditions. An abrupt 
reduction in leakage under the highest speed condition in Figure 5-12 occurs when the rate 
of dissipation is just sufficient to boil the fluid within the body of the seal. The stiffness 
curves for the lowest three speeds in Figure 5-12 exhibit strong positive stiffness due to the 
dominance of the inlet loss effect. The stiffness curve for the highest speed in that plot 
undergoes an abrupt change at the outset of boiling and so has a different character than the 

others. 

Figures 5-13 and 5-14 give leakage information for a constant speed of 30,000 RPM 
for various film thicknesses and degrees of subcooling. Figure 5-13 shows that the 
leakage rates for a highly subcooled initial state are much higher than those for moderate 
subcooling. In Figure 5-14, the curve of the highest degree of subcooling shows strong 
positive stiffness. In that curve, the flow was pure liquid throughout the body of the seal 
so the effect of the inlet losses were dominant. Note that the opening forces for the large 
film thicknesses is much lower than their moderately subcooled case counteiparts. Among 
the moderate subcooling cases, the relative order of opening force in terms of degree of 
subcooling shifts with film thickness for reasons similar to those given for the curve shifts 
in speed in previous examples. Most importantly, this plot shows that small changes in the 
degree of subcooling in the upstream reservoir can drastically affect the seal behavior 
changing the stiffness from positive to neutral or negative at a given mean film thickness. 
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Figure 5-13: Effects of film thickness and subcooling on 
6 leakage rate. 

Sealed fluid is Oxygen. 

Rotation Speed, <0 = 30,000 RPM 
Reservoir Pressure, P<» = 4.3 MPa 
Back Pressure, Pb = 0.0 MPa 

Coning Angle, p = 0.0 radians 




Figure 5-14: Effects of film thickness and subcooling on 

opening force. 

Sealed fluid is Oxygen. 

Rotation Speed, (0 = 30,000 RPM 
Reservoir Pressure, Poo = 4.3 MPa 
Back Pressure, Pb = 0.0 MPa 

Coning Angle, p = 0.0 radians 
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Figures 5-15 and 5-16 give leakage rates and opening forces for a seal with various 
angles of coning. Note that coning not only gives a small decrease in leakage, but also 
gives the seal a substantial positive stiffness over the entire range of film thicknesses. 

Figures 5-17 and 5-18 give leakage and stiffness information with water as the sealed 
fluid under various degrees of subcooling. The character of the water seal behavior is 
similar to that of the cryogenic oxygen seal. Note that the force-displacement curves for the 
two highest degrees of initial subcooling show substantial positive stiffness. For those 
curves, under the given operating conditions, the flow was all-liquid throughout the seal 
for every film thickness shown. 

Consideration of some trial parametric studies indicate the following salient features 
of turbulent face seal behavior: 

• Vapor production through pressure drop and viscous dissipation within the seal 
promotes choking and reduces the leakage rate. 

• Subcooling the incoming fluid partially negates the effects of dissipation and 
pressure drop on leakage. 

• Centrifugal inertia effects tend to retard flow in outside seals and reduce leakage. In 
inside seals, the opposite is true. 

• All-liquid choked flow may occur at low rotational speeds and high degrees of inlet 
subcooling. The leakage rates associated with this type of choked flow may be 
significandy larger than those seen in two-phase flows. 

• The interplay of the effects of dissipation, subcooling, and rotational speed is 
complicated and may give rise to peculiar effects. The character of the seal opening 
force and stiffness may change radically for small changes in the operating 

parameters. 
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Unexpected positive stiffness (i.e., a force-displacement relation that gives a 
restoring force) is found in all-liquid flow situations. Decreasing the mean film 
thickness reduces the inlet losses and so increases the opening force. 

Coning has a beneficial effect on seal leakage and promotes positive seal stiffness. 
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Figure 5-15: 


Effects of film thickness and coning on 
leakage. 

Sealed fluid is Oxygen. 


Rotation Speed, to = 30,000 RPM 
Reservoir Temperature, T«, = 150.0 K 
Reservoir Pressure, P« =4.3 MPa 
Back Pressure, Pb = 0.0 MPa 
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Figure 5-16: Effects of film thickness and coning on 
6 opening force. 

Rotation Speed, co = 30,000 RPM 
Reservoir Temperature, Too = 150.0 K 
Reservoir Pressure, Poo = 4.3 MPa 
Back Pressure, Pb = 0.0 MPa 

Coning Angle, p = 0.0 Radians 
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Figure 5-17: Effects of film thickness and subcooling on 

leakage. 

Sealed fluid is Water 
Rotation Speed, (0 = 2,500 RPM 
Reservoir Pressure, Poo = 3.16 MPa 
Back Pressure, Pb = 0.100 MPa 

Coning Angle, (J = 0.0 Radians 
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Figure 5-18: Effects of film thickness and subcooling on 
opening force. 

Rotation Speed, (0 = 2,500 RPM 
Reservoir Pressure, Poo = 3.16 MPa 
Back Pressure, Pb = 0.100 MPa 

Coning Anble, (3 = 0.0 Radians 
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CHAPTER 6 

LEAKAGE IN TWO-PHASE ANNULAR SEALS 
6.1 Introduction 

The quasi-steady analysis of the turbulent annular seal roughly parallels that of the 
radial face seal. Upon applying the integral approach to the basic equations, equations 
similar to those governing the face seal flow are obtained. The "jump" conditions relating 
the reservoir thermodynamic state to the inlet state are identical with those of the previous 
problem and the general solution method proceeds as before. Moreover, leakage rate 
solutions to annular seals will be seen to have the same general character as those for radial 

face seals. 

This seal will be idealized as having polished working surfaces which may be 
considered hydraulically smooth. Further, the seal shaft will be taken to be centered in the 
bearing surface so that the geometry under study is axisymmetric. For reference, a 
schematic diagram of a typical annular seal is given in Figure 6-1. Small displacements of 
the shaft from this centered position are not expected to change the leakage greatly so that 
the leakage characteristics of annular seals under real operating conditions may be estimated 
reasonably well by the present model. Of course, no forces lateral to the long axis of the 
shaft will be developed if the shaft is in its centered position. 
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Figure 6-1: Annular seal geometry (not to scale). 


It will be seen that the choking, a phenomenon which limits the leakage rate, is 
sensitive to the rate of internal heat generation through viscous shear and to the inlet 
thermodynamic state. The previous studies of Childs, et al. [127, 128, 129) do not 
consider the possibilities of two-phase flow, nor do they account for internal heat 
generation. This analysis then represents a more complete accounting of the important 
effects influencing the leakage rate than was previously done. 

6-2 Specialized Equations for Liquid, Two-Phase, and Vapor Flow 

Dropping the terms dependent on the circumferential location from the basic equations 
given in section 4-2 gives the following: 


Continuity: 

2jtRhG = m (6-1) 
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Axial Momentum: 

0 di + ^ + ^2 = 0 

dz dz h 


Adiabatic Energy Equation: 


di 

dz 



Rco . 

- Gh T rel R 


( 6 - 2 ) 


( 6 - 3 ) 


It should be noted that the continuity equation requires no specialization; it holds in its 
present form regardless of the phase of the flowing fluid. Also, the circumferential 
momentum equation is trivially satisfied as before. The mean circumferential velocity is 
assumed to be fully developed throughout the interior of the seal with the value of w = 
(0R/2 since the film thickness, and hence the circumferential flow development region, are 
very small compared to the other characteristic dimensions of the seal. 
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Liquid Flow 


As before, the liquid density is taken as a constant with a value corresponding to the 
density of saturated liquid at the same temperature. Equation (6-1) then implies that the 
axial velocity is constant. This affords further simplification of the axial momentum 
equation (6-2) which becomes: 

dp . 2 (6-4) 

dz h rz 

and of the energy equation (4-2-3) which becomes: 

di Rg> , (6-5) 

to “ Gh^ lR 

These equations are always well behaved, giving a monotonically decreasing pressure 
and a monotonically increasing enthalpy. Since the liquid compressibility is neglected, the 
sonic speed for the liquid is infinite. Normal choking, where the fluid velocity reaches the 
acoustic wave speed in the medium, then can not be predicted from the equations above. It 
should be noted that the acoustic wave speeds for most liquids are so high that choking in a 
mode as described above is very seldom experienced. 

Liquid - Vapor Flow 

The two-phase flow through the seal is assumed to be homogeneous so that the fluid 
properties are constant through the film. The axial derivatives of quality and enthalpy are 
expressed in terms of the axial derivative of pressure as was done earlier. The axial 
velocity in equation (6-2) and (6-3) may be expressed in terms of the mass velocity, G. 
Using this and the relations among the saturation properties, these equations may be solved 
for the axial derivative of pressure explicitly, assuming G is known: 
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( 6 - 6 ) 


dp = 

dz 


1_ 

h (1 - <t>) ! 


2 1 + 


*fg 

G 2 yv t g, 


lt n + 



X 




where, 


<t> = - 

T y 


a. 


(1 + G^/Bpl*) ifg 

G w fg 


and, the axial gradient of mixture quality may be simply expressed in terms of the pressure 
gradient as: 

<& _ 2T rz (1 +G 2 av/aplx) dp (67) 

^ = ' hG 2 v fg ' G 2 v fg ** 


Since the temperature is only a function of pressure at saturation, the following is true: 


di_ 

dz 


(di f ^dp - dX 
\dP K dPjdz + lf «dz 


( 6 - 8 ) 


where the pressure gradient is known from (6-6) and the quality gradient is known from 
(6-7). Note that equations (6-6) and (6-8) have a singularity at conditions where <|> equals 
unity. It is seen that the equations governing two-phase flow in this case are very similar to 
those of the face seal case and that the nature of their singularities are also similar. Thus, 
the two-phase, turbulent annular seal can exhibit normal and all-liquid choked behavior 

much like that of the radial face seal. 
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AH - Vapor Flow 


If one assumes an ideal gas with constant specific heats, the continuity, axial 
momentum, energy, and state equations may be combined in a simple way to give the 
gradients of the pressure and enthalpy explicitly. The basic equations adapted for vapor 

flow are then: 
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and. 


di 

dz 
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1 + (y- 1) M 



A/i 2 d P 
+ TMv^ 


( 6 - 10 ) 


This is a well known one-dimensional gas flow problem with heating and friction. Here, 
M is the Mach number. 

6.3 Numerical Examples and Discussion 

A nominal design for the interstage seal of the Space Shuttle Main Engine High 
Pressure Oxidizer Turbopump is taken to be: 

Shaft Radius, R = 0.0325 m 
Seal Length, L = 0.0260 m 
Radial Clearance, h = 1.74 x 10 -4 m 


and the range of operating conditions is taken as; 

Rotation Speed, (0 = 20,000 to 30,000 RPM 
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Reservoir Temperature, Too — 116 to 143 K 
Reservoir Pressure, P» = 2.75 to 4.83 MPa 
Back Pressure, Pb = 0 to 0. 1 MPa 

These values will be used unless stated otherwise. The flow inlet is taken to be square- 
edged with a nominal inlet loss coefficient of 0.5. 

Figure 6-2 shows temperature, pressure, and quality profiles for a seal under 
unchoked conditions in the two-phase regime. Although the back pressure in this instance 
is out of the design range, it is useful to present this case for comparison with others. 
Here, the reservoir fluid is sub-cooled by 1 K. Note that boiling had been initiated at the 
inlet so that the quality there was non-zero. The profiles in this case are smooth with 
moderate gradients throughout the seal. The leakage rate in this case was 0.575 kg/s. 

Figure 6-3 shows the profiles for choked conditions in the two-phase regime. The 
reservoir fluid is sub-cooled by 5 K and boiling in this case occurs well inside the seal. 
Discontinuities in the slopes of the pressure and temperature profiles occur at the location 
where boiling begins. Unlike the case before it, the pressure gradient at the exit has an 
extremely large negative value. The leakage rate in this case was 0.686 Kg/s. 

Figure 6-4 shows the profiles for the anomalous case of all-liquid choked flow. The 
reservoir fluid is sub-cooled by 24 K. The profiles shown are nearly straight lines and the 
leakage rate is 1.34 Kg/s, almost twice that of the other cases. Although the fluid exit 
velocity is much lower than the sonic speed of the liquid, the mass flow is at its maximum 
value and is independent of back pressure, provided it is below a critical value. 

The following three numerical examples are studies of the effects of various 
parameters on the leakage rate. The geometry and operating conditions will be taken to be 
the same as in the case presented previously in Figure 6-3 unless stated otherwise. 
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Figure 6-5 shows .he dependence of leakage rate on die length to radius ratio on the 
shaft rotational speed. As either rotation speed or path length increase, more heat is 
dumped into an element of fluid as it flows through the seal. This increased heating brings 
about increased vapor generation. The mixture density is reduced, and die axial velocity is 
increased. Compressibility effects, however, limit die flow so .ha, the net mass leakage 

rate decreases. 
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Fieure 6-2: Unchoked two-phase flow through the seal. 

Shaft Radius, R = 0.0325 m 
Seal Length, L = 0.026 m 
Radial Clearance, h = 1.74 x 10 -4 m 
Shaft Rotation Speed, <a = 30,000 RPM 
Reservoir Temperature, Too = 139.0 K 
Reservoir Pressure, Poo = 2.79 MPa 
Back Pressure, Pb = 1-80 MPa 
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Figure 6-3: Choked two-phase flow through the seal. 

Shaft Radius, R = 0.0325 m 
Seal Length, L = 0.026 m 
Radial Clearance, h = 1.74 x 10 -4 m 


Shaft Rotation Speed, (0 = 30,000 RPM 
Reservoir Temperature, Too = 135.0 K 
Reservoir Pressure, Poo = 2.79 MPa 
Back Pressure, Pb = 0 MPa 
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Fieure 6-4: Choked all-liquid flow through the seal. 

Shaft Radius, R = 0.0325 m 
Seal Length, L = 0.026 m 
Radial Clearance, h=1.74xl0' 4 m 
Shaft Rotation Speed, co = 30,000 RPM 
Reservoir Temperature, Too = 1 16.0 K 
Reservoir Pressure, Poo = 2.79 MPa 
Back Pressure, Pb = 0 MPa 
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Figure 6-5: Effect of UR ratio and rotation speed on leakage. 

Shaft Radius, R = 0.0325 m 
Radial Clearance, h = 1.74 x 10" 4 m 
Shaft Rotation Speed, 0) = 30,000 RPM 
Reservoir Temperature, T« = 135.0 K 
Reservoir Pressure, Poo = 2.79 MPa 
Back Pressure, Pb = 0 MPa 


Figure 6-6 shows the dependence of leakage rate on the radial clearance and on the 
shaft rotational speed. Small film thicknesses and high rotational speeds promote high 
rates of heat generation by viscous dissipation. As in the previous example, the leakage 
rates associated with high rates of heat generation are small. 

Figure 6-7 shows the dependence of leakage rate on the degree of sub-cooling in the 
upstream reservoir and on the shaft rotational speed. Sub-cooling at the upstream reservoir 
partially negates the effect of the viscous heat generation in the film and so brings about 

higher leakage rates. 
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Fieure 6-6: Effect of film thickness and rotation speed on leakage. 

Shaft Radius, R = 0.0325 m 
Seal Length, L = 0.026 m 

Shaft Rotation Speed, 0) = 30,000 RPM 
Reservoir Temperature, T«, = 135.0 K 
Reservoir Pressure, P«» = 2.79 MPa 
Back Pressure, Pb = 0 MPa 

In summary, consideration of some trial parametric studies show that the leakage 
characteristics of turbulent annular seals are similar to those of radial face seals discussed 
earlier. Vapor production through pressure drop or viscous heat generation in the seal is an 
important mechanism for limiting the leakage rate. Small radial clearances, high rotation 
speeds, and extended seal lengths all aid in reducing leakage. Sub-cooling in the upstream 
reservoir inhibits vapor production and causes large leakage rates. High degrees of sub- 
cooling in the upstream reservoir may bring about all-liquid flow, even at very low back 
pressures. Leakage rates for all-liquid flows can be much higher than two-phase or vapor 

flows and should be avoided. 
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Fieure 6-7: Effect of sub-cooling and rotation speed on leakag - 
Shaft Radius, R = 0.0325 m 
Seal Length, L = 0.026 m 
Radial Clearance, h = 1.74 x 10- 4 m 
Shaft Rotation Speed, ca = 30,000 RPM 
Reservoir Pressure, P«, = 2.79 MPa 
Back Pressure, Pb = 0 MPa 
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CHAPTER 7 


A COMBINED LAMINAR AND TURBULENT COMPUTER 
CODE FOR FACE SEAL PERFORMANCE PREDICTION 
AND DESIGN 


7 . 1 Introduction 


A computer code which combines the discrete boiling laminar model developed in 
Chapter 3 and the adiabatic turbulent model developed in Chapter 5 is presented here. 
Using this code, we can calculate the performance of a face seal under both low and high 
leakage conditions. Once we determined how this seal performs under the two limiting 
conditions, we are able to interpolate the results and predict how the seal will behave in the 
intermediate region, where both models fail. We performed a parametric study using the 
computer code to investigate the effects of various operation parameters on face seal 
performance. The parameters being investigated are the level of subcooling, coning of seal 
faces, rotational speed, conductivity of the seal materials and the width of seal faces. 
Because our analysis is based on idealized models, the numerical results may not be exactly 
the same as those for an actual seal, where a complete analysis must involve consideration 
of the complex geometry of the seal, its surroundings and other complicated effects such as 
misalignments of axes and surface roughness. The results presented here are intended to 
be viewed as a rough guide and are used to indicate the general trends of how seal 
performance is affected by varying these parameters. The combined laminar and turbulent 
code is documented in Appendix D. 


7.2 Instability of Face Seals 

Liquid near saturation condition is notorious in its difficulty to be sealed 
successfully. Seal instability resulting from phase change has been identified as the cause 
of numerous failures as for example in the Forties Main Oil Line pump seals [84], 
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Harrison and Watkins [84] observed .hat the Forties Main Oil Line pump seal failures 
usually occurred after the pumps had been operating at above normal temperatures and 
when the water content in the crude oil was high. Under these operating conditions, fluid 
vapor was present in a large portion of the seal as vaporization started early on in the 

leakage path through the seal. 

In previous studies using single component fluids as the working fluids [87, 88, 
89, 90], the authors show that it is possible for a two phase face seal to have two 
equilibrium operating points, one stable and the other unstable. The stable operating point 
is at a higher f.lm thickness and its boiling location is closer to the seal exit than the 
unstable point. It is shown that if phase change started very close to the inlet, a seal will 
exhibit negative stiffness and be unstable. This occurs when the sealed liquid is near 
saturation. Although the fluid in the case of the Forties Main OB Line [84] was a mixture 
many components, it is conceivable that tire mechanisms that caused the instability of of the 
pump seals were quite similar as the operating temperature was close to the saturation 

temperature of water. 

There are also reports of failures of reactor coolant pump seals during nuclear 
power station blackouts which were brought about by similar mechanisms [82, 102], 
Normally, these seals operate with well-cooled water under high pressure. During a stanon 
blackout, external cooling to the coolant system is lost and sometimes there is also a drop in 
the system pressure. The coolant water can quickly reach near saturation condition and 
cause the seals to become unstable. The reactor coolant pump seals may "pop open" under 
these conditions, resulting in uncontrollable and excessive leakage which in severe cases 

can lead to a loss of coolant accident. 

Our results show that, for a parallel seal, there is an absolute minimum amount of 
subcooling which must be provided to maintain stable operation. A seal will become 
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unconditionally unstable if the temperature of the sealed liquid rises above this limit. It is 
shown that stability can be increased by positive seal coning, where the gap between the 
seal plates converges in the direction of the flow, but at the expense of higher leakage. 
However, depending on the balance, "popping open" of seals may occur at a temperature 
lower than the absolute stability limit. It is because at temperatures below the absolute 
stability limit, any increase in the fluid temperature will lead to an increase in the seal 
opening force. If the closing force was not high enough to balance the increased opening 
force, the seal will pop open. The temperature at which this would occur is determined by 

the balance ratio of the seal. 

7.3 Discussion of Sample Calculations and Seal Stability 

This study was performed using water as the working fluid. Sample calculations 
including stiffness curves for parallel and coned face seals are presented to show the effects 
of the various operation and seal parameters. 


7.3.1 Effects of Subcooling 


The parameters chosen for the sample calculations are: 


p 0 =2000kPa 
pj = lOlkPa 
k = 50W/(m K) 
coning slope = O.Om/m 


r 0 = 0.0428625m (1 ll/16in.) 
ri = 0.0365125m (1 7/16in.) 
speed = 4000rpm 


To illustrate how subcooling affects the operation of a face seal by changing the 
boiling location, we look at the pressure profile of the fluid along its leakage path in a 
typical parallel outside face seal. Figure 7-1 shows plots of pressure profiles along the 
leakage path for an all liquid, a two-phase and an all vapor seal under low leakage 
conditions. When subcooling is high, the fluid remains liquid from inlet to exit. In most 
seals where the widths of the faces are small compared to their radii, the pressure drop is 
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approximately linear from inlet to exit (see curve A in Figure 7-1). For a parallel face, low 
leakage seal, where inlet loss and inertial effects are unimportant, the opening force is a 
function of the pressure differential between the inlet and exit only. Thus, it is quite simple 
to obtain a good estimate of the opening force of a parallel all liquid face seal. This all 
liquid opening force is always less than the opening force generated by the same seal under 
two-phase or all vapor operation. With a sufficient decrease in subcooling (i.e. the sealed 
liquid approaches saturation), the liquid starts to boil at the exit and the boiling location 
gradually moves towards the inlet as subcooling is reduced. Pressure drop in the vapor 
region is much steeper than in the liquid region of the seal, as it is shown in Figure 7- 1 (see 
curve B). As the boiling location move from the exit towards the inlet, the seal opening 
force initially increases rapidly, reaches a maximum and then decreases eventually to the all 
vapor seal opening force, which is always higher than the all liquid opening force. Curve 
C in Figure 7-1 shows the pressure profile of an all vapor seal. If the leakage rate is high 
and the flow is turbulent, the inlet pressure loss may be significant and if choking occurs, 
there is a sudden drop in the pressure at the outlet. Hence, depending on the conditions, 
the opening force of a turbulent seal may be higher or lower that the equivalent low leakage 
laminar seal. 

It should be noted that once the boiling location has moved past the point where the 
opening force is maximum, the seal is unstable. This instability can be explained by 
looking at the axial stiffness of the seal. When the boiling location moved past the point at 
which the opening force is maximum, any decrease in film thickness due to external 
disturbances will cause the boiling location to move further towards the inlet because of the 
increased heat generation. This reduces the opening force. The seal plates will move 
closer together and the seal will eventually collapse. If the seal is perturbed to move apart, 
the opening force will increase. Depending on the level of subcooling and the balance 
ratio, the seal will either blow open or reach the stable operating point which is at a higher 
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Him thickness. Therefore, for a parallel face seal, we can see that there is a minimum 
amount of subcooling we must provide to maintain its stability, regardless of the seal 

balance. 

In actual operations, the closing force of a seal changes very little with the seal film 
thickness. As subcooling varies, the film thickness adjusts itself such that the opening 
force is in equilibrium with the closing force. As an example, Figure 7-2 shows a stiffness 
plot (plot of opening force versus film thickness) of a parallel seal for different degrees of 
subcooling. In this example, the saturation temperature of the sealed liquid is 485.6K. If 
we assume a balance such that the seal closing force is constant at 2500N, the seal will 
operate stably at a very small film when the reservoir liquid is highly subcooled (see Figure 
7-2). For example, the operating film thickness for 390K bulk liquid temperature of is less 
than 10' 7 m. Surface contacts are likely to occur at such a low film thickness. As 
subcooling is reduced, the stable operating film thickness increases and as a result, the 
leakage rate also increases. In Figure 7-2 we show, when the bulk fluid temperature nses 
above 460K (more than 25K below the saturation temperature of the fluid in the reservoir), 
the is no stable operating film thickness predicted by the laminar seal model which gives an 
opening force of 2500N, and the seal pops open. Figure 7-3 shows the plot of leakage rate 
versus film thickness and Figure 7-4 shows the plot of seal stiffness coefficients versus 
film thickness. Notice that when the bulk fluid temperature is above 480K, the laminar seal 
model predicted that the seal would be unconditionally unstable. 

Since the quasi-isothermal discrete boiling laminar model is valid only for small 
leakage rate seals, the sample calculation results given become more and more inaccurate as 
the film thickness gets higher and higher. Using the computer code developed based on the 
adiabatic turbulent model, we may determine the seal behavior at large film thicknesses and 
high leakage rate. The adiabatic model assumes turbulent flow and negligible heat 
conduction into the seal plates when compared to heat convection within the fluid. 
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Convection terms in the fluid and inlet pressure loss are considered and boiling is allowed 
be continuous in this model. Full account is taken for all liquid, two-phase and vapor 
choking. Figure 7-2 shows the sample results of the parallel seal using both the quasi- 
isothermal laminar model and the adiabatic turbulent model. When subcooling is high, the 
the seal opening force given by the adiabatic turbulent model decreases with increasing film 
thickness. This is due to the increased inlet pressure loss as the radial fluid velocity 
increases with film thickness. As subcooling reduces, the opening force increases 
dramatically as the flow starts choking at the exit. The trends of how opening force varies 
with film thickness predicted by both models are very similar. When the laminar model 
predicts that the opening force drops as film thickness increases, the turbulent model 
predicts a similar behavior. When the laminar model predicts that the seal opening force 
increases as film thickness increases, the turbulent model predicts that as film thickness 
increases even further, the seal opening force will increase even more. The intermediate 
region, where both heat conduction into the seal plates and heat convection within the fluid 
are important, is bounded quite well by the two limiting models. This close agreement of 
the predicted behavior of seals using both models gives us confidence in applying the 
laminar model to predict seal failures. 

Figure 7-5 shows the absolute minimum reservoir temperature for stable operation 
of a parallel face seal as determined using the laminar model. When rotational speed is 
zero, a low leakage seal with negligible inlet loss is always neutrally stable because the 
opening force is constant with film thickness. The opening force of a two phase stationary 
parallel face seal is determined by the boiling location and is a function of the reservoir fluid 
temperature alone. The reservoir temperature at which the opening force of this seal is 
maximum is the highest temperature the seal can maintain stable operation. We can reason 
this by looking at a seal running at an infinitesimally slow speed. If the seal temperature is 
at the stability limit, any decrease in the film thickness will generate more heat, thereby 
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moving the boiling location closer to the inlet and lowering the seal opening force. Since 
the seal closing force remains constant, the seal plates are pushed together until the seal 
collapses. This temperature limit calculated from an infinitesimally slowly rotating seal can 
be considered as the absolute stability limit. It represents the best case and is by no means 
a conservative estimate. Any seal running at an appreciable speed will have a different 
temperature limit which will always be at a lower temperature than the one denved using 
the stationary seal. Also, in order to maintain successful seal operation, the seal 
temperature often can not exceed a temperature limit which is much lower than the stability 
limit. This operating temperature limit is dependent on the seal balance ratio. If the seal 
balance ratio chosen is low, the increase in seal opening force due to the reduction of 
subcooling can pop open a seal and result in a seal failure. Using the models and computer 
code presented here, one can calculate this limit for a seal given the balance ratio. A line 
representing such a limit is shown in Figure 7-5 for an arbitrary seal balance ratio. 

7.3.2 Effect of Coning 

Figures 7-6 through 7-11 show results of sample calculations for seals with coning 
(converging in the direction of flow) using the same operation parameters and seal 
dimensions as in the last section. Two cases of seal coning are presented. Figures 7-6 
through 7-8 show the results of a seal with a minute coning slope (10- 5 m/m) and Figure 7- 
9 through 11 a larger coning slope (1.5xl0-3m/m). For the seal with the smaller coning 
slope, the effects of coning are not very noticeable at the larger film thicknesses (see Figure 
7-6) and the opening force and leakage rate characteristics are quite similar to that of the 
parallel seal (see Figures 7-6 and 7-7). When the film thickness is small, coning effectively 
restricts the flow and thereby increases the seal opening force and stability (see Figures 7-6 
and 7-8). As the film thickness gets very small, there is little pressure drop in the liquid 
region of the seal. When the film thickness decreases further, unlike a parallel face seal, 
the boiling location moves back closer to the seal exit and the opening force asymptotically 
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approaches the maximum opening force the seal can produce due to hydrostatic pressure. 
This maximum opening force can be calculated by taking the fluid pressure inside the seal 
be uniformly at the reservoir fluid pressure. The asymptotic limit can be explained by 
looking at the limiting case that, when the plates of a coned seal just touch each other at the 
exit, there is no leakage and if the centrifugal effects are not important, the fluid inside the 
seal gap will be uniformly at the reservoir pressure. For high speed seals, the maximum 

opening force will be lower. 


For the seal with the smaller coning slope, the film thickness at the seal exit is only 
a minuscule 6 x 10 * 0 , smaller than at the inlet but the effects of coning are already 
significant. Any wear of the seal plates at the outer radius resulted from surface contacts 
will produce a much higher coning slope. Figures 7-9 through 7-1 1 show the results of a 
seal with a more realistic coning slope. The difference between the inlet and exit film 
thicknesses is about 10 =m. If we assume the same seal balance as in the previous 
example, the stable operating film thickness is much higher for the coned face seal than for 
the parallel face seal. At 450K bulk fluid temperature, the equilibrium film thicknesses are 
about 4xl0' 7 m for the parallel seal (see Figure 7-2) and 3x10-%. at the inlet for the coned 
seal (see Figure 7-9) and the predicted leakage rates are lxlO ^gfs and 3xl0'kg/s 
respectively. The high leakage rate of the coned seal, however, invalidated the 
assumptions made in the laminar seal model. Using the results from the adiabatic turbulent 
seal model, the equilibrium film thickness predicted is a lower 1.7x10-%. and the leakage 
rate is about 6xl0- 7 kg/s (the results of the leakage rate calculations using the adiabatic 
turbulent model are not shown). Since the leakage rate is very high in this case, the 
predictions using the turbulent model should be more accurate. A leakage rate this high is 
unacceptable for most sealing requirements and therefore, in order to maintain an acceptable 
leakage rate, we should choose a higher balance ratio for a coned seal than we would for a 
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parallel seal. At higher temperatures, the stable operating film thicknesses are so high that 
the seal can be considered popped open. 

It is interesting to note that, unlike the parallel face seal, the unstable operating 
region of a coned seal is either very small or non-existing. Since the seal opening force 
drops when the temperature exceeds the value at which the opening force is maximum, the 
seal stiffness characteristics at a temperature very close to saturation can be quite similar to 
that of a much lower temperature (see Figure 7-9). It is possible that a seal with a certain 
balance ratio can operate satisfactorily when subcooling is high or very low but not in 
between. However, once a seal is popped open at a lower temperature, it is not very likely 
that the seal will regain stable operation by itself. Figures 7-10 and 7-1 1 show the leakage 
rate and stiffness coefficients as predicted using the laminar seal model; because of the 
higher operating film thicknesses, the stiffness coefficients of the coned seal are much 
lower than those of the parallel seal (note the change of scale of the vertical axis). 

7.3.3 Effect of Speed 

Figures 7-12 through 7-17 show the stiffness curves for seals running at speeds 
from 0 rpm to 10000 rpm and at constant bulk fluid temperatures. The basic operation 
parameters and seal dimensions used in the calculation were the same as in the previous 
two sections. Figures 7-12 through 7-14 show the results of a parallel seal operating at 
bulk liquid temperatures of 400K, 450K and 470K respectively. Figure 7-15 through 7-17 
show, under the same operating conditions, the results of a coned seal with a converging 
coning slope of 1.5xlO- 3 m/m. From Figures 7-12 through 7-14 we can see, for the 
parallel seal with a given seal balance, as the speed increases, the stable operating point 
shifts towards higher film thicknesses and as a result, the leakage rate also increases. This 
is because when the seal speed increases, the rate of heat generation, and therefore the seal 
temperature, increases. This moves the boiling location closer towards the inlet and 
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increases the seal opening force when the seal is operating in the stable region. Thus for a 
given balance, the stable equilibrium operating point must move to a higher film thickness 
when the seal speed is increased so that it can maintain the same opening force (and also a 
similar boiling interface location). For the turbulent region, the heat generation rate is 
relatively small and the effect of speed is not very significant 

It can be seen from the figures that, even at 10000 rpm, the centrifugal inertia 
effects are not very significant. However, if either the seal speed gets much higher or the 
pressure differential between the inlet and exit becomes much lower, the centrifugal effects 

can become important. 

For the coned seal, where heat generation cannot be very high because of the seal 
geometry, the effects of speed is not significant in both the laminar and the turbulent 

regions (see Figures 7-15 through 7-17). 

7,3.4 Effects of the Thermal Conductivity of Seal Materials 

Figures 7-18 through 7-25 show the stiffness curves of seals made with materials 
of different thermal conductivities. In our studies presented in previous sections, the 
average conductivity of the seal materials used is 50W/m-K. To illustrate the effects of 
varying thermal conductivity, analyses are performed using the same operating parameters 
and dimensions for seals made with lower and higher thermal conductivity materials. 
Figures 7-18 and 7-19 show the stiffness plots of two parallel seals for different seal 
speeds with average seal material conductivities of 15W/m-K and lOOW/m-K respectively. 
Figures 7-20 and 7-21 show the stiffness plots of the same seals for different bulk 
temperatures. It can be seen from these curves (and from Figures 7-13 and 7-2), for seals 
made of lower conductivity materials, the stable equilibrium operating film thicknesses are 
higher. This is because for the seals with lower thermal conductivities, under the same 
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operation conditions for the same film thicknesses, the seal temperatures is higher than the 
ones with higher conductivities and thus they have higher opening forces. 

Figures 7-22 through 7-25 show stiffness curves for different temperatures and seal 
speeds for coned seals (coning slope 1.5xl0- 3 m/m) with average thermal conductivities of 
15W/m-K and lOOW/m-K. Since the heat generation rates are relatively small, the effects 
of seal material conductivity is not very significant. 

For the high leakage turbulent region, since our model is adiabatic, the values of 
seal material thermal conductivity have no effect on the results. 

7.3.5 Effects of the Width of Seal Faces 

Figures 7-26 through 7-29 illustrate the effects of varying the widths of the seal 
faces. Figure 7-26 shows the stiffness curve for a parallel seal with a face width half of 
that of the seal used for the studies in the previous sections and Figure 7-27 shows the 
curve for a seal with a face width twice as wide. From these curves we can see that the 
opening forces is almost direct proportional to the seal face areas. When compared to the 
results of the original seal shown in Figure 7-2, the stability characteristic for all three seals 
are very similar to each other (i.e. the shape of the stiffness curves are very similar) and the 
temperatures at which the seal become unstable are almost identical. 

Figure 7-28 and 7-29 show the stiffness plots for coned seals with face widths half 
and double the width of the seal used in the studies in the previous sections. Again, when 
compared to the results of the original coned seal (see Figure 7-9), the effects of the width 
of seal faces on the stability characteristics are shown not to be significant. 
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7.4 A Summary of Our Results and a Discussion of the Implications to 
Face Seal Designs 


Using the computer code base on the quasi-isothermal low leakage model and the 
adiabatic high leakage model, we can determine the stability criteria and stiffness 
coefficients of face seals. Our results show that these important seal operation 
characteristics predicted by the two extreme models tend to overlap each other in a 
reasonable manner. This gives us confidence that the simplified quasi-isothermal model is 
capable of being a useful tool for face seal design and can be used to form the basis of a 
face seal design methodology. 

Subcooling is a critical factor in determining seal stability. The simplified model 
allows us to predict very quickly how the seal opening force and stiffness characteristics 
vary with subcooling. As subcooling of the sealed fluid is reduced, the seal opening force 
increases rapidly. It is shown in out results that, even when the seal stiffness coefficient is 
positive (i.e. the seal is operating stably), given a sufficient reduction in the subcooling of 
the sealed fluid, the seal may still pop open if the balance force is exceeded. Using the 
model and the computer code presented here, one can establish the temperature limit for 
successful seal operation for different inlet pressures given a seal balance ratio, such as the 

example shown in Figure 7-5. 

The effects of seal coning has also been discussed. In actual seal operations, any 
surface contact caused by wobbling of the seal plates tends to wear the plates at the outer 
radius. A parallel face seal can become coned very quickly once it enters into service. It is 
shown in our sample calculation results that a small amount of seal coning can dramatically 
alter the operation characteristics of a face seal, in particular its stiffness characteristics. A 
coned seal is in general more stable than a parallel face seal (compare Figure 7-2, 7-6 and 
7-9). Although coning tends to stabilize a seal, given the same balance ratio, a coned seal 
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operates at a higher film thickness and thus has a much higher leakage rate than a parallel 
seal. It should be realized that when determining the seal balance ratio, a higher value is 
needed for a coned seal than for an otherwise equivalent parallel seal. Besides the wear of 
seal plates, during actual operations, effective coning can also be caused by temperature 
and pressure distortion. An equivalent amount of coning can be added to account for these 

effects. 

The effects of speed, thermal conductivity of seal materials and the width of the seal 
faces are also discussed. It is found that, under normal operating conditions (i.e. high 
pressure differential between the inlet and exit and seal speed much lower than the critical 
speed), the centrifugal inertia effects are not very significant. As seal speed increases, 
given the seal balance, the stable equilibrium operating point moves towards higher film 
thicknesses and the leakage increases. It is also shown that, for the same balance, the 
equilibrium operating film thickness of a seal made with lower thermal conductivity 
materials is higher than that of a seal made with higher conductivity materials. Both these 
effects can be explained by looking at the seal temperature. When the seal is operating in 
the stable region, an increase in seal temperature due to any reason, such as an increase in 
the bulk fluid temperature, an increase in seal speed or because the seal is made fan lower 
thermal conductivity materials, will increase the seal opening force. It is found that, the 
width of seal faces has relatively little effects on the temperature limit for successful and 
stable seal operations and die seal opening force is almost directly proportional to the area 

of the seal faces. 

To ensure a long and reliable face seal operation, a seal designer must know the 
level of subcooling of the sealed fluid, anticipate how much it will vary and take into 
account the effects of coning resulting from wear and pressure and thermal distortions, and 

then the balance ratio can be chosen accordingly. 
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The computer code presented here currently incorporates the physical and 
thermodynamic properties of water, nitrogen, oxygen and hydrogen. Other fluids can be 
easily added into the code with only minor modifications. Also, the solution to the heat 
transfer problem is based on an idealized model in which the seal area is surrounded by two 
semi-infinite solids. If higher accuracy of the seal temperature profile is needed, a 
geometry specific heat transfer calculation, such as a finite element analysis, can be 
incorporate into the present code. 

Although no experimental results are presented here, published field observations 
and some unpublished experimental data tend to corroborate our results. 
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Figure 7-1: 


Plot of Pressure Profiles of a Typical Face Seal Under 
All Liquid, Two-Phase and All Vapour Operations 
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Figure 7-2: 


Stiffness Curves of a Parallel Face Seal 
For Different Bulk Fluid Temperatures 
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a Parallel Face Seal for Different Bulk Fluid Temperatures 
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Figure 7-8: 


Plots of Stiffness Coefficients versus Film Thickness 
of a Face Seal with a Small Coning Slope for 
Different Bulk Fluid Temperatures 
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Figure 7-11: Plot of Stiffness Coefficients versus Film Thickness 
of a Coned Face Seal with a Typical Coning Slope 
for Different Bulk Fluid Temperatures 


168 



Opening Force (N) 


4500 r 


4000 l 


3500 


3000 


2500 


2000 


1500 


1000 


500 


ADIABATIC TURBULENT MODEL 
ISOTHERMAL LAMINAR MODEL 


/ 


/ 

/ S~ 

///-" 
/ / / ■ — 

Zv/Z 

ZZZZ- 

//"Z'/'Z 

'///Z V// 


2000 rpm 
3000 rpm 

4000 rpm 
5000 rpm 
€000 rpm 
7000 rpm 
8000 rpm 
9000 rpm 
10000 rpm 


//yyVy 

/ / s ' '// / / / / 

4«y 
/ / //,■ / / // 



\ 




- 

Bulk Tsmparatura 

400.00 

K > 


Conductivity 

50.00 

W/»-K 


Inlat Radius 

0.042862 

m 


Out la t Radius 

0.036513 

m 


Coning Slops 

0 . 00 a +00 

m/m 


Inlat Prassura 

2000.00 

kPa 


Outlat Prassura 

101.00 

kPa 


Saturation Tamp. 

485.58 

K 


10000 rpm 
0 rpm 


1.00+08 


1 . 00+07 


1 . 00+06 


1 . 00+05 1 . 00+04 1 . 00+03 

Film Thickness at Inlet (m) 


Figure 7-12: Stiffness Curves of a Parallel Face Seal for Different 
Seal Speeds at 400K Bulk Fluid Temperature 
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Figure 7-13: 


Stiffness Curves of a Parallel Face Seal for Different 
Seal Speeds at 4S0K Bulk Fluid Temperatures 
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Figure 7-15: Stifness Curves of a Coned Face Seal for Different 
Seal Speeds at 400K Bulk Fluid Temperature 
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Fieure 7-16* Stiffness Curves of a Coned Face Seal for Different 
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Fieure 7-17: Stiffness Curves of a Coned Face Seal for Different 
Seal Speeds at 470K Bulk Temperature 


174 


o 



Opening Force (N) 


4500 r 


4000 


3500 


3000 


2500 


2000 \ 


ADIABATIC TURBULENT MODEL 
ISOTHERMAL LAMINAR MODEL 



0 rpm 

1000 rpm 
2000 rpm 
_ 3000 rpm 
• ^ 4000 rpm 

5000 rpm 
, _ 6000 rpm 

7000 rpm 
8000 rpm 
9000 rpm 
^ 10000 rpm 


7500 1 


10000 rpm 
0 rpm 


000 “ Bulk Temperature 

Conduct i vit y 
Inlet Radiue 
Outlet Radiue 
SOQ Coning Slope 

Inlet Free rare 
Outlet Freeeare 
Saturation Temp . 

1.So^08 1.00+47 1.00+46 


450.00 

K 

15.00 

W/m-1 

0.042862 

m 

0.036513 

m 

0.00*4-00 

m/m 

2000.00 

kFa 

101.00 

kPa 

485.58 

K 


> -i — J M 

1.00+45 1.00+44 1.00+43 

Film Thickness at Inlet (m) 


Figure 7-18: Stiffness Curves of a Parallel Face Seal with a 

Low Thermal Conductivity for Different Seal Speeds 
at 450K Bulk Fluid Temperature 
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Figure 7-20: 


Stiffness Curves of a Parallel Face Seal with a Low 
Thermal Conductivity for Different Bulk Fluid Temperatures 
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Figure 7-22: 
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Figure 7-23: 


Stiffness Curves of a Coned Face Seal with a 
High Thermal Conductivity for Different Seal Speeds 
at 450K Bulk Fluid Temperature 
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Figure 7-26: 


iffness Curves of a Parallel Face Seal with a Narrow 
ice Width for Different Bulk Fluid Temperatures 
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Figure 7-27: 


Stiffness Curves of a Parallel Face Seal with a Wide 
Face Width for Different Bulk Fluid Tempreatures 
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Figure 7-29: 


Stiffness Curves of a Coned Face Seal with a Wide 
Face Width at Different Bulk Fluid Temperatures 
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Chapter 8 


A Summary of Detailed Work and Key to Publications 

We have discussed some of the distinctive behavior characteristics of two-phase 
seals, particularly their axial stability. While two-phase seals probably exhibit instability to 
disturbances of other degrees of freedom such as wobble, etc., under certain conditions, 
such analyses are too complex to be treated at present. Since an all liquid seal (with parallel 
faces) has a neutral axial stiffness curve, and is stabilized axially by convergent coning, 
other degrees of freedom stability analyses are necessary. However, the axial stability 
behavior of the two-phase seal is always a consideration no matter how well the seal is 
aligned and regardless of the speed. Hence, we might think of the axial stability as the 
primary design consideration for two-phase seals and indeed the stability behavior under 
sub-cooling variations probably overshadows other concerns. The main thrust of this 
work has been the dynamic analysis of axial motion of two-phase face seals, principally the 
determination of axial stiffness, and the steady behavior of two-phase annular seals. 

A chart is shown below which identifies the problems and we have considered and 
indicates location of the detailed analyses and computer codes in the listed publications. 


187 



Type of Program Description 
Seal 


Application 


Reference 
for details and 
computer codes 


Face Sea] 
High Leakage 
Laminar 

Adiabatic Laminar Steady, 

Stability, 
Dynamic Response, 
Limit Cycles 

High leakage 
laminar seal, 
not encountered 
under normal 
operating conditions. 

Papers: 2, 4 
Thesis 1 

Face Seal 
Low Leakage 

Quasi- Isothermal Laminar: 

Steady, 
Stability, 
Elevated temper- 
ature calculated 
with heat transfer 
model 

Good for most 
face seals under 
normal conditions 
with low leakage 

Papers: 2, 6 
Thesis 1 
Chapter 7, 
this report 

Face Seal 
Low and Moderate 
Leakage 

Continuous boiling with 
with variable temperature. 
Laminar film coefficients 
considered in boiling region, 
homogeneous T-P model. 
Convection included 
in energy equation: Steady, 
Stability, 
Transient 

Refinement of 
above model. 

Good for virtually 
all low and moderate 
leakage face seals. 

Paper: 8 
Thesis 3 

Face Seals 
Turbulent 
High Leakage 

Adiabatic turbulent Steady, 

face seals including Stability 

inertia, heat generation 
and entrance losses 

High leakage seals, 
special applications 
such as cryogenic 
pumps: LOX-GOX 

Papers: 3, 5 
Thesis 2 
Chapter 7, 
this report 

Annular Seals 
Turbulent 
High Leakage 

Adiabatic turbulent Steady, 

annular seals Stability 

including inertia, heat 
generation and entrance 
losses 

Cryogenic LOX- 
GOX Turbo Pumps 
for Rocket Engines 

Papers: 3, 7 
Thesis 2 
Chapter 7, 
this report 

General 

Design 

Codes 

A simplified combined computer 
program for laminar and turbulent 
face seals including effects listed 
in items 2 and 4 above 


Chapter 7, 
this report 
(disk inside 
back cover) 


188 




List of published references cited in above Table: 

1. N.S. Winowich, M.J. Birchak, W.C. Kennedy, and W.F. Hughes, Phase Change 
in Liquid Face Seals," ASME Paper No. 77-LUB-12, Trans. ASME, Journal of 
Lubrication Technology, Vol. 100, No. 1, January 1978, p. 74. 

2. N.H. Chao and W.F. Hughes, "Phase Change in Liquid Face Seals D-Isothermal 
and Adiabatic Bounds with Real Fluids,' Transactions of the ASME, Journal of 
Lubrication Technology, Vol. 102, No. 3, July 1980, p. 350. 

3. R.M. Beeler and W.F. Hughes, "Turbulent Two-Phase Flow in Ring and Face 
Seals," Proceedings of Ninth International Conference on Fluid Sealing, BHRA 
Fluid Engineering, April 1-3, 1981, pp. 185-202. 

4 R M Beeler and W.F. Hughes, "Dynamics of Two-Phase Face Seals," Trans. 
A.S.L.E., Vol. 27, No. 2, April 1984, p. 146. 

5 PA Beatty and W.F. Hughes, "Turbulent Two-Phase Flow in Face Shaft Seals," , 
Journal of Tribology, Trans, of the A.S.M.E., Vol. 109, No. 1, January 1987, pp. 
91-99. 

6 P. Basu, R.M. Beeler, and W.F. Hughes, "Centrifugal Inertia Effects in Two-Phase 
Face Seal Films" Trans. A.S.L.E., Vol. 30, No. 2, April 1987, pp. 177-186. 

7. P.A. Beatty and W.F. Hughes, "Turbulent Two-Phase Flow in Annular Seals" . 
Trans. A.S.L.E., Vol. 30, No. 1, January 1987, p. 11. 

8. P. Basu and W.F. Hughes, "Thermal Instability in Two-Phase Face Seals" , 
Proceedings of 11th International Conference on Fluid Sealing, BHRA, Cannes, 
France, 1987, pp. 423-441. 

Ph.D. Theses: 

1. Beeler, R.M., 'Effects of Two-Phase Flow in Shaft Seals,’ Ph.D. Thesis, Carnegie 
Mellon University, May 1985, Pittsburgh, Pennsylvania. 

2. Beatty, P.A., ’Inertial Effects in Turbulent Fluid Seals,’ Ph.D. Thesis, Carnegie 
Mellon University, May 1986, Pittsburgh, Pennsylvania. 

3. Basu, P., ’Thermal Effects in Two-Phase Flow Through Face Seals,’ Ph.D. Thesis, 
Carnegie Mellon University, May 1988, Pittsburgh, Pennsylvania. 
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10. APPENDICES 


Appendix A 

Derivation of the Energy 
Equation 


Consider radial viscous flow with properties uniform across the film. The tem- 
perature is assumed to vary only with r (but small variations across the film 
account for gradients in that direction and consequent conduction into the seal 
rings). The energy equation is written below, 



Ti + * + * v 


Jj 


where k is the thermal conductivity of the fluid and $ is the dissipation function 
given to a high degree of accuracy by 



Integrating the energy equation across the film 


r-i* - 

£*(£♦£*£)* 


Multiplying the r-momentum equation (2.1) by u and integrating across the film, 





r 
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Combining with the energy equation and neglecting conduction along the film 

m di dv h f h titr* f h /0u’\ J , t h , 

— — — = fiv— + p I dz + / p I — I dz + k—dz 

2ir r dr dz 0 Jo r Jo \ dz ) Jo d 3 z 

Since 1 / = 0 at z = 0,h, the first term of the right-hand side goes to aero. The 
term 


/>£ 


dz 


integrates to 


k 


dT 

dz 


z=h 



~9 


where q is the conduction heal flux into the seal rings from the fluid. 

Using the expression for t/ and w from equations [2.4] and [2.3], the following 
form of the integrated energy equation is obtained, 

m di 

2tt r dr ^ h 


3u? 2 m p 7 r 7 h*u? A 
20* + 700ji 
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Appendix B 

The Steady State Influence 
Coefficient Matrix 


In this appendix, the relationship between the heat-conduction profile and the 
temperature profile is derived. The result is a coefficient matrix, [C], which 
relates the temperature and heat conduction vectors, {T} and {?}, whose n 
elements correspond to finite difference points that have been placed along the 
fluid film. The elements of this matrix involve the complete elliptic integrals of 
the first and second kind, K and E , respectively. These are defined as follows : 



where r is referred to as the modulus. 

It is required to find the temperature T at at a point on the surface of the seal 
ring faces due to a heat source at another position. The background temperature 
is Toe. The radial position r, at which temperature is to be determined, corre- 
sponds to the mid-width of the » th element of the fluid film. The heat source 
q t is due to uniform heat generation over the f** element. The heat flux leaves 
the fluid film and enters the rings normal to the faces. The heat flux is assumed 
constant across each element and is thus discontinuous at the edge between two 
adjacent elements. The inner and outer radii of the i ,h element are designated 
as r,_ and r t+ respectively. 

The temperature T at due to a point heat source q, at rj in the solid is 


T(r-„r t ) 



202 



which iE easily derived by considering the heat flux across a sphere of radius 
— rj| centered around a point heat source of intensity q t at f t . The temperature 
T tt is now given by integrating the above Green’s function over the entire i th 
element as follows: 


T. t 



rq t d<f>dr 
Airk |r, - r t | 


+ r oc 


Carrying out this integration, according to the formulas provided by Abramowitz 
and Stegun [106], and Byrd and Friedman [107], the following expressions are 
obtained, 



Equation (B.l), (B.2), or (B.3) is used depending on whether the temperature 
in question is at a radial position outside, inside or equal to the corresponding 
position of the heat source, respectively. Since, in general, all of the elements 
making up the fluid film will be releasing or absorbing heat, it is summed over 
all the heat sources to get the total temperature elevation at a point. 


T. 


n 


£r. t + :r oe 

1=1 


The temperature T, and the heal fluxes q t are grouped into column vectors {T} 
and {9}, respectively, which are then related by the square coefficient matrix [C], 

{ T } = [C]{q} 4- {7U 


where {Too} is the column vector, each element of which is equal to the back- 
ground temperature T, 
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Appendix C 

Derivation of the 
Time-Dependent Influence 
Coefficient Matrix 


Referring to the Figure (C.l), the temperature at the mid-point of the element 
‘A:’ due to heat generation over the element */’ over the time interval (0,/) is given 
by the following Equation [108], 


Tki(1) = 


2nqiriAria 


i: 


1 

*■*+*■? 1 

I f rfcr< 

[*a(f - t )] 3/J ^ 

4q(< — t) 

* [2a(t - r) 


dr (C.l) 


where, A r/ U the width of the ‘/th’ element and 7o is the modelled Bessel 
function. Here a time invariant concentrated line ring source is assumed to exist 
at the radius Vj\ 


C.l Off-Diagonal Terms (k ^ /) 

Evaluate the integral J\. 

Assume id = i — r. Hence dtj = —dr. For different ranges of different 
approximations are to be used. 
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c.1.1 Case 1 

P° r *“'* > ^ or <d < = (**)« ( sa y), the following series representation of 

Bessel function /„ can be made according to the Equation (9.8.2) in [106], 

r = 6j + h.o.i (C.2) 

where b. , = 0.39894228 and higher order terms (h.o.t) can be neglected. 

Using the Equation (C.2), the integral J, can be written as 


»i /“* 1 [ (r*- r,)*l 

31.499 o^vT/o eT7, [ 4^~J^ 

p f ( r * ~ ri)* j 

31.499 a^rOri [ 4 of, J < C * 3 ) 


where E y is the ‘Exponential Integral Function 1 the 
which is given by the Equation (5.1.53) in [106]. 
Substituating Equation (C.3) in (C.l), 


series representation of 


r "(°,M = 0.199471 qt 

= ^w(0,fj)9; ; 0< t d <(<c)*/ (C.4) 

where C k i is the (*,/)th element of the time-dependent influence coefficient 
matrix. 


C.l. 2 Case 2 

For (f e )u < t d < oo t the modefied Bessel function I„ can be expressed by the 
following senes expression according to the Equation (9.8.1) in [106] with the 
values of the coefficients <to through a 4 given, 



Using the Equation (C.5) in evaluating the integral Jj, one obtains 
= 8tt 3 / ? o y~ » ~ 4 - | [3*5449 /jo(y) + 0.56889 ai-Rju/n(y) + 


0.16182 aj/?w/u(y) + 0.04603 a 3 /?*/ 13 (y) + 0.0131 m/&/n(y)]* 

l*(y)]£ (0.6) 


1 


8ttVJ q 


\/ r 2 + r ? 


where y(y) is the function inside the parentheses and the limits of integration 
are given by 


yi = 


yj = 


d ± r l 

4q<j 

rl ± if 

4q/» 


The other functional definitions, used in the Equation (C.6) are 

= + *) 

/io(y) = er/( v /y) 

/n(y) = 1.32934 er/( v ^)-(y + 1.5) v ^e-*' 

/u(y) = 11.63173 er/( v /y)-(y 8 + 3.5y l + 8.75y + 13.125) v /ye-*' 
fn{y) = 287.8853 er/( y'y) - (y 8 + 5.5y 4 + 24.75y* + 86.625y* + 216.5625y + 
324.8438) v /ye~ v 

/i 4 (y) = 14034.407 er/( v /y) - (y T + 7.5y # + 48.75y # + 268.125y 4 + 1206.563y 8 + 
4222.969y J + 10557.422y + 15836.13) N /y«~'' 

Substituating Equation (C.6) in (C.l) with the limits of integration, 

■ 


= C w (f„<j) 9 , ; (< e ) w < U < oo 

Again qi is assumed to be constant over the time interval / j and f 2 . 


(C.7) 
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C.2 Diagonal Terms (k = /) 

Considering an heat source extended over the ‘Arth’ ring, the temperature at the 
mid-point of the ring ‘A:’ due to self heat generation over the period (0, <) is given 
by 


Tkh{0,i) = 

r f r r > i 


n: 


2 rrq k 1 

ri+r> 1 

T 

r k r 

pC 8 [»a(f - t)]* /s eIP 

1 

iU 

0 ^ 

1 

*0 

2 a (t — t) 


rdrdr 


(C.8) 


Again let t d = t-r. Using this substituation in Equation (C.8), there obtains, 


Tkk(0,i) 


2irq k a f 1 1 [ T n 

Jo Jr.. " P 



rdrdij 


(C.9) 


C.2.1 Case 1 


For 0 < I 4 < (/ e )*k , using the series representation (C.2) of the modeled Bessel 
function, the Equation (C.9) can be written as 


Tkt,(0,1) = 


qkb 1 M r*+ 

>/8rrr fc k Jo tj «/r t _ ^ 


4 at< 


y/rdrdtj 


(C.10) 


The function y/r inside the spatial integral doe6 not vary that much over the 
interval (r* + ,r*_ ) and it can be assumed as y/r* over the same interval. However, 

the first function exp ~ j varies dramatically over the region specially when 

u -* 0. 

Hence the Equation (C.10) becomes, 


rr , n ^ r 1 r k + ~ r ) 

r " (0 ' ,) 7mJ.uL'* r r~Mr\ 


Ul Ul^ 


*1 

k 


I ft * V vg " r k — ^ 'XUVg J 

(™) 4 ^ "'(&)]* 

- . n v j ^ / j \ 


= Ckk 9 * ; 0 < < (/ f )*ir 


(C.ll) 


C.2. 2 Case 2 

For (/ f )kk < < 00 , the same Equation (C.7) as in A: / / can be used. 
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cn:m?2 r r 



*c 


Figure C.l: Discretization of Annular Seal Interface 
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Appendix D 


The Combined Laminar and Turbulent 
Low and High Leakage Operations 


Computer Code for 


D.l Input File Description 


The default name of the input file is 'leakin'. If a file named 'leakin' is found in 
e cunent working directory, the program will assume this file to be the intended input file 
^ WUi ^ ^ it without asking the user any question. IfluTa m 

n °' eiUS ' " ,h ' *• - « prompted for the name of Ute input Z 


The format of an input file is tabulated below: 


Line 

1 

2 

5 — 

Variable 

outfnm 

fluid 

Format 

al5 

a8 

Description 

fjame ol the output data file, maximum 15 characters 

However, for machines running MS-DOS this is 

touted by the MS-DOS file natong convention 
( character file name + 3 character extension) 

iName of the working fluid, ‘lbe current fluids 
incorporated in this program are Vater^mWen' 
hydrogen’ and 'oxygen'. nitrogen , 

4 

5 

6 

pinlet 

pback 

rinlet 

f20.0 

f20.0 

hf20.0 " 

~ at the ■»!« (l.e.reservoir pressure! 

assure at the seal exit (..<=, back pressure). ' 

Radius of tne seai inlet hor an outside seal, this is the" 
outer radius; for an inside seal, this is the inner radius. 

V 

n 

routit 

i20.0 

“Kaoius of the seal outlet. For an outside seal, this is — 

radius 11 ^ radlUS; f ° r 8,1 lnSidc Sea1, 11115 is * e outer 
fm] 

/ 

ft 

cone 

f20.0 ' 

a^H^° ning S, °P?‘. A P° sitlve value means that the seal' 
seal increasin g racial distance from the 

[m/m] 

u 

conauc 

f20.0 

[W/m-Kl COnaUCtlVlty ° f thC 8681 face maten als. 
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Line Variable 




Format 





next 

ncurve 

lines 


tinfs 

or 

rpms 


f20.0 


Description 


o decide whether the analysis will be performed for 
both low leakage and high leakage operations or only 
for the low leakage operation. 

If a ’Y' or y is entered, the analysis will be 
performed for both flow regimes. If any other 
character is entered (an TV or ’n' is recommended) the 
analysis will only be performed for the low leakage 
laminar operation. 

Note that a combined analysis for both flow regimes 
typically takes 40 tunes (depending on the number of 
different film thicknesses calculated in the high leakage 
analysis 1 ) longer to run than an analysis for low 

leakage operation only. 

To choose whether the analysis will be repeated for 
different seal speeds or for different bulk (reservoir) 
fluid temperatures. 

If an 'S' or ’s' is entered, the analysis will be repeated 
for different seal speeds. 

If a T' or 't' is entered, the analysis will be repeated 
for different bulk temperatures. 


the analysis is going to be repeated for different seal 
speeds (choose = ’S'), the bulk temperature tinf is 
entered here. 

If the analysis is going to be repeated for different bulk 
temperatures (choose = T'), die seal speed rpm is 
entered here. 

[K] if choose - 'S'; [rpm] if choose = T 

Number of times the analysis will be performed for 
different seal speeds or bulk temperatures. This is the 
same as the number of stiffness and leakage curves to 
be calculated. 

The limit on the number of curves is 50. 


The ncurve different seal speeds (if choose = ’S') or 
bulk fluid temperatures ( if choose = T') for which 
analyses are performed are entered here. 

Each seal speed or bulk fluid temperature is entered on 
a separated line. 

[rpm] if choose = 'S'; [K] if choose = T* 


Section D.3 describes how to vary the range and the number of film thicknesses. 
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A sample input file is given below. Columns 21 to 80 of each input line are 
reserved for comments. 


seal. out 

water 

2000. 

101 . 

0.0428625 

0.0365125 

1.5e-03 

50.0 

y 

t 

4000.0 
10 

360.0 

390.0 

420.0 

450.0 

460.0 

470.0 

475.0 

480.0 
482.5 

485.0 


output file name (max 15 characters, less for DOS) 

fluid (either: water, nitrogen, oxygen or hydrogen) 

pinlet (kPa) 

pback (kPa) 

rinlet (m) 

rout It (m) 

cone (m/m) 

conduc (W/m-K) 

both, ’n' if laminar analysis only; 'y' for both 

choose, 's' for speeds, 't* for bulk temperatures 

if 's', bulk temp, (degrees K); seal KPM if 't' 

nuiber of curves (max 50) 

tinf (1) 

tinf (2) 

tinf (3) 

tinf (4) 

tinf (5) 

tinf (6) 

tinf (7) 

tinf (8) 

tinf (9) 

tinf (10) 


The ranges of allowable temperatures and pressures for the various fluids currently 
incoiporated in this program are listed below: 


Water 

305.00 K 

< 


4.718 kPa 

< 

Nitrogen 

63.15 K 

< 


12.54 kPa 

< 

Oxygen 

80.00 K 

< 


30.09 kPa 

< 

Hydrogen - 

13.80 K 

< 


7.042 kPa 

< 


tinf 

< 

647.29 K 

pinlet, pback 

< 

22.089 MPa 

tinf 

< 

126.20 K 

pinlet, pback 

< 

3.400 MPa 

tinf 

< 

154.58 K 

pinlet, pback 

< 

5.043 MPa 

tinf 

< 

32.94 K 

pinlet, pback 

< 

1.284 MPa 
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D.2 Output File Description 

The contents of output file is self descriptive. An output file begins with a heading 
showing the input seal operating parameters. It is followed by tabulated results of the seal 
analyses for each seal speed or bulk fluid temperature. The tabulated results are all given in 
S.I. units. If an analysis is performed for both low and high leakage operations, the results 
for each seal speed or bulk temperature will be tabulated in two tables, the first table 
contains the results from the quasi-iso thermal laminar model and the second table contains 
the results from the adiabatic turbulent model. If the analysis is performed only for low 
leakage operations, the results will be tabulated in only one table. 

The output data contains : 

For each seal speed or bulk fluid temperature, a sub-heading is printed to identify either 
the speed or the bulk temperature: 

rpm: Seal speed, [rpm]; or 

tinf: bulk fluid temperature of the reservoir, [K]. 

For the laminar low leakage model, the tabulated results are: 
hrin: The inlet film thickness, [m]. 

xboil: The non-dimensional boiling location. The non-dimension radial 

position, x, is defined as fraction of the seal width from the seal exit: 

r - r exit 

X = 

r inlet" r exit 

pboil: The fluid pressure at the boiling locations, [Pa], 

tseal : The temperature of the seal at the boiling location, [K]. 

leakage: The leakage rate [in kg/s], where positive values indicate that the leakage 
path is radially outwards and negative values indicate that the leakage 
path is radially inwards, 
load: The seal opening force, [N]. 

• Note that if the xboil and pboil columns are blank, the seal is either all-liquid or 
all-vapor. For an all-liquid seal, the temperature tseal is the temperature 
evaluated at the seal exit. For an all-vapor seal, since we assumed that the 
viscous dissipation is negligible, tseal will be the same as the bulk fluid 
temperature. 
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• If the pboil and leakage columns are blank, the seal is non-boiling and is 
operating at above critical speed In this case, the non-dimensional radial location 
of the stationary liquid front is given as xboil and the leakage rate is zero. 

For the turbulent high leakage model, the tabulated results are: 


hrin 
Reyc: 
Reyr: 
qexit. 
leakage 


load: 


The inlet film thickness, {m]. 

The maximum Reynold's number in the circumferential direction. 

The maximum Reynold's number in the radial direction. 

The quality of the fluid exiting the seal. 

Tfie leakage rate [in kg/s], where positive values indicate that the leakage 

path is radially outwards and negative values indicate that the leakage 
path is radially inwards. 

The seal opening force, [NJ. 

D.3 Preset Program Parameters 

Th, are several parameters preset (hard coded) in the computer code. Some of these 
ptcameters are problem specific and may be adjusted to suit the specific needs of individual 

cod T- t, ■ eyarc lmnal ‘ Ie ‘ J " I ' med,a,dy “P°" Program execution. The location of the 
«.e whtch tmualizes these parameters is in the main program following the variable 

declanmons. A bnef description of each of these parameters are tabulated below: 



closs 


epsi 


real*8 


real *8 


0.95 


“-Pressure toss coefficient. The inlet pressure 
toss is given by the product of the inlet loss 
coefficient and the dynamic pressure. The preset 
value is for inlets with square edges 

una contraction factor. The factor by which the" 
finite difference gnd contracts along the leakage 
path. Since pressure drop is usually much 
steeper near the exit, we can obtain better 
resolution using a variable sized grid with grid 
points closer to each other near the exit 
The iteration e rror bound for the exit p n»cc„„.' 


213 



Preset parameters for the number and range of film thicknesses to be calculated: 


Variable 

Type 

Preset 

Value 

Description 

npntst 


■ 

Number of points (different film thicknesses) on 
the turbulent curves. The points are chosen in 
such a way that they are equally spaced between 
hinupt and hinlot in log scale. 

hinupt 

real* 8 


The upper limit of the inlet film thickness of the 
turbulent curves. 

[m] 

hinlot 

real*8 

10-5 

The lower limit of the inlet film thickness of the 
turbulent curves. 

[m] 

npntsl 

integer*4 

II 

Number of points (different film thicknesses) on 
the laminar curves. The points are chosen in 
such a way that they are equally spaced between 
hinupl and hinlol in log scale. 

hinupl 

real*8 

10-5 

The upper limit of the inlet film thickness of the 
laminar curves. 

[m] 

hinlol 

real* 8 

io - 8 

The lower limit of the inlet film thickness of the 
laminar curves. 

[m] 


In the original turbulent code, the user is required to input floor, an estimate of the 
minimum allowable mass flow rate. This value is used as the first mass flow rate estimate 
when solving the boundary value problem using the shooting method. If this estimate is 
too large and gives a choked flow in the first iteration of the shooting method, the program 
will abort (return without a solution) with an error message. If the estimate is too low, the 
program may not converge. 

Fortunately, the range of values of floor for which the program will converge is quite 
large. Since the mass leakage rate is directly proportional to the film thickness, in the 
combined program, the value of floor is initialized in the turbulent analysis code 
subroutine face. The current preset values are, where h mcan is the mean film thickness: 


floor 


10 2 kg/s, if h mcM1 > lO^m 
10 4 kg/s, if h mean < lO^m 
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Note that the optimal choices for the values of floor are problem specific and the preset 
values may not be suitable for all cases. If the turbulent code does not run properly, the 
user should try adjusting the preset values of floor. 
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D.4 Source Code of the Combined Laminar and Turbulent Program 



Steady State Analysis for a Two-Phase Face Seal Operating in 
Low Leakage Laminar and the High Leakage Turbulent Regimes 


This program analyzes the steady state performance of a 
two-phase face seal operating under both low leakage 
conditions using the quasi-isothermal laminar model arKl 
under high leakage operations using the adiabatic turbulent 

model . 


lit S I t 

S.I. units are used in all calculations. With the exception 
of input pressures being in kPa, all input values are in 
S.I. units. All output values are in S.I. units unless 
they are stated otherwise. 


Acceptable Ranges of Temperature and Pressure Input Data 


Water 

- 305.00 K 

< 


4.718 kPa 

< 

Nitrogen 

- 63.15 K 

< 


12.54 kPa 

< 

Oxygen 

- 80.00 K 

< 


30.09 kPa 

< 

Hydrogen 

- 13.80 K 

< 


7.042 kPa 

< 


tinf 

< 

647.29 K 

pinlet, pback 

< 

22.089 MPa 

tinf 

< 

126.20 K 

pinlet, pback 

< 

3.400 MPa 

tinf 

< 

154.58 K 

pinlet, pback 

< 

5.043 MPa 

tinf 

< 

32.94 K 

pinlet, pback 

< 

1.284 MPa 


Please direct any questions or comments to: 

Professor William F. Hughes or Stephen Lau 
Department of Mechanical Engineering 
Carnegie Mellon University 
Pittsburgh, Pennsylvania 15213 





c+ c 

c . C 

c Main program. 

c 

program seal 
c 
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c + 


+c 


implicit logical ( a - z ) 
implicit undefined ( a - z ) 

real*8 r inlet, routlt, cone, conduc 

real*8 pinlet, pback, rpm, omega, rolinf, tint 

real *8 closs, epsi, error t 

real *8 hinupt, hinlot, himpl# honlol 

real* 8 hinlet, pboil, rboil 

real *8 tseal, visgas, visliq, rolliq, rgas 


cc 

cc 


integer nspace 
integer npntst, npntsl 


common /geomty/ r inlet, 
conmon /operat/ pinlet, 
common /turbin/ nspace, 
common /curves/ npntst, 
common /iter at/ hinlet, 
conmon /fluid / tseal, 


routlt, cone, conduc 
pback, rpm, omega, rolinf, tinf 
closs, epsi, error _ 

hinupt, hinlot, npntsl, hintpl, hinlol 

pboil, rboil 

visgas, visliq, rolliq, rgas 


cc 


real*8 rpms(50), tinfs(50) 
real *4 tarray (2) , dtime 
integer i, ncurve 
logical exstf 
character fluid* 8, respon*l 
2 both*l, dummy* 60 


inpfnm*15, outfnm*15, choose*!, 


c — 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


The built-in turbulent program parameters are 
Seme of these parameters are problem specific 
according to individual problems. 


initialized below, 
and may be adjusted 


NSPACE - 
CLOSS - 

EPSI 


ERROR - 


ize of the finite difference grid. 

-ilet loss coefficient. The inlet pressure loss 
is given by the product of the inlet loss 
coefficient and the dynamic pressure, 
rid contraction factor. The factor of which the 
grid size contracts along the leakage path. 

Since pressure drop is usually much steeper near 
the exit, we can obtain better resolution using 
a variable sized grid with grid points closer to 

each other near the exit. 

-it- oration error bound for the exit pressure. 


cere is one other preset parameter FWOR ; an ft urate of the 
■ninimum allowable leakage rate. In the original _ f^lf t 
Droaram written by Paul Beatty, the user is required to inpot 
FLOOR. This value is used as the starting est urate m solving 
the boundary value problem using the shooting method. If the 
SfufS^OOR is too high and gives a choked flow in the first 
iteration of the shooting method, the turbulent routine will 
abort (return without a solution) with an error message. If 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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FLOOR is too low, the program may not converge. . 

Fortunately, the range of values of FLOOR for which the program 
wUl^on^^Tis^ite wide. Since the leakage rate is directly 

related tothe fita thickness, the value of FtfCR is irntial- 
ized in the turbulent analysis SUBROUTINE FACE. The current 
preset^values aie, «here HMEfiN is the mean fita thickness: 


c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c 


FLOOR = 1.0d-02 [kg/s] ; if HMEAN > 1.0d-04 

= 1.0d-04 [kg/s] ; if HMEAN <= 1.0d-04 


Note that the optimal choices for the values of FLOOR 
problem specific and these preset values may not be suitable for 
all cases. If the turbulent part of the code does not run 
properly, the user should try adjusting the preset values o 

FLOOR. 


C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c 


[m] 

[m] 


nspace = 50 
closs = 0.5d0 
epsi = 0 . 95d0 
error = 1 . Od-3 


c FLOOR is initialized in the turbulent analysis code 
C FACE. _ 
c ~~ 


c 

SUBROUTINE C 

c 

c 


c 


c — 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c~ 


The~upper and lower limits of film thicknesses in the study and 
the number of points are decided below. They may be changed to 
suit specific needs. 


NPNTST - Number of points (different film thicknesses) on 
the turbulent curves. The points are chosen in 
such a way that they are equally spaced between 
HINUPT and HINLOT in log scale. 

HINUPT - The upper limit of the inlet film thickness of 
the turbulent curves (in meters) . 

HINLOT - The lower limit of the inlet film thickness of 
the turbulent curves (in meters) . 

NPNTSL - Number of points (different film thicknesses) on 
the laminar curves. The points are chosen in 
such a way that they are equally spaced between 
HINUPL and HINLOL in log scale. 

HINUPL - The upper limit of the inlet film thickness of 
the laminar curves (in meters) . 

HINLOL - The upper limit of the inlet film thickness of 
the laminar curves (in meters) . 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

— c 


npntst =21 
hinupt = 1.0d-03 
hinlot = 1.0d-05 
npntsl = 31 
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hinupl = 1.0d-05 
hinlol = 1.0d-08 


c 7 c 

c Open the input data file. Default name of the input file is c 
c ' seal .in'. c 
c 


inquire (f ile= ' seal .in', exist =exstf) 
if (exstf) then 

open ( 1 1 , f ile= ’ seal . in status= ' old ' ) 
else 

write (*,'(" File "seal . in" does not exist . . . " ) ' ) 
write (*, ' (' 1 Please enter the input file name' ’, 

2 " (max. 15 characters) : ")') 

read(*, ' (al5) ') inpfnm 
open ( 11 , file=inpfnm, status= ' old' ) 
end if 

c c 

c Read the inputs . . . 

c 

c 

c For some reason, Microsoft FORTRAN compiler (at least V3.31) 
c requires the read statement to read in the comment field along 

c with the relevent input data on each input line or the program 

c bombs out with an I/O error. To get round this, the comment 

c field of each line is read in as a character string DUNMY. 

c c 

read (11, 1 (al5,5x,a60) ') outfnm, dummy 
read (11, ' (a8,12x,a60) ') fluid, dummy 
read(ll, ' (bn,f20.0,a60) ') pinlet, dummy 
read ( 11 , ' (bn, f20.0,a60) ') pback, dummy 
read(ll, ' (bn,f20.0,a60) ') rinlet, dunmy 
read(ll, ' (bn, f20.0,a60) ') routlt, dummy 
read(ll, ' (bn, f20.0,a60) ’) cone, dummy 
read(ll, ' (bn, f20.0,a60) ') conduc, dummy 
read(ll, ' (bn,al, 19x,a60) ') both, dummy 
read(ll, ' (bn,al, 19x,a60) ') choose, dummy 

cc read (11, ' (al5) * ) outfnm 

cc read ( 11 , * (a 8 ) ') fluid 

cc read ( 11 , ' (f 20 . 0 ) ') pinlet 

cc read ( 1 1 , ' (f 20 . 0 ) ') pback 

cc read ( 11 , ' (f 20 . 0 ) ') rinlet 

cc read ( 11 , ' (f 20 . 0 ) ' ) routlt 

cc read ( 11 , ' (f 20 . 0 ) ') cone 

cc read ( 11 , ' (f 20 . 0 ) ') conduc 

cc read ( 11 , ' (al) ') both 

cc read (11, ' (al) 1 ) choose 

if ( (choose. eq. ’S') .or. (choose. eq. * s’) ) then 
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oooooo on 



cc 

cc 


cc 

10 


cc 

cc 


cc 

20 


read (11» ' (bn, f20.0,a60) ') tinf, dumtiy 
read (11, ' (f20.0) ' ) tinf 

read (11, ' (bn, il0,l0x,a60) ') ncurve, dummy 
read (11, 1 (ilO) ') ncurve 
do 10, i = 1, ncurve 

read (11, ' (bn, f20.0,a60) ’) rpnsd), dummy 
read(ll, ’ (f20.0) ') rpms(i) 

else 0 ?? 1 Uchoose . eq. ' T ' ) .or. (choose. eq. 't') ) then 
read (11, ' (bn,f20.0,a60) ' ) rpm, dummy 

read(ll,'(f 20.0)') rpm 

read (11, ' (bn, il0,10x, a60) ') ncurve, dunmy 
read(ll, ' (ilO) ') ncurve 
do 20# i = 1/ ncurve 

read (11, ' (bn, f20.0,a60) ') tinfsd), dummy 

read (11, ' (f20.0) ') tinfsd) 

continue 

el «rite(*, ’ (" Invalid input for CHOOSE, ' 

• 1 program terminates « • • ) / 

stop 
end if 


close (11) 


c 

c 
c 

c information 
c 





from the fluid properties data 


call rsatdt (fluid) 

c Open the output file ... _ 

c Microsoft FORTRAN doesn't know about thefile^sta^ £ 


cc 


inquire (f ile=outfnm, exist=exstf ) 

if ^ited^M" File ” ’ ' , al5, ' ' " already exists, type", 
’ it i » i t y» • » » to overwrite it, ' ' ) ' ) outfnm 

write (*, ' (48x, ' ' else to exit : ' ' ) * ) 
read(*, ' (al) ') respon _ ^ 

if ( (respon. eq. "i') .or. (respon. eq. y )) then 
open (12, file^outfnm) 
open (12, file=outfnm, status- unknown ) 

else 
stop 
end if 

else 
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open (12, file=outfnm, status='new' ) 
end if 

c c 

c Write the input seal parameters as the heading of the output file c 
c so as to identify the output file. c 

c c 

if ( (choose. eq. 'S' ) .or. (choose. eq. ’s') ) then 

write (12, 501) pinlet, pback, rinlet, routlt, cone, conduc, 

2 tinf, fluid 

501 format (lx, 'Inlet Pressure = ',0pfl0.4,4x, ' kPa' / 

2 lx, 'Exit Pressure = ',0pfl0.4,4x, ’ kPa' / 

3 lx, 'Inlet Radius = ',lpel4.4, ' m' / 

4 lx, 'Exit Radius = ',lpel4.4, ' m' / 

5 lx, 'Coning Slope = ',lpel4.2, ' m/m' / 

6 lx, 'Conductivity = ',0pfl0.4,4x, ' W/m-K' / 

7 lx, 'Bulk Temperature = ',0pf8.2, 6x, ' K' / 

8 lx, 'Fluid = ' , 3x, a8 /) 

else 

write (12,502) pinlet, pback, rinlet, routlt, cone, conduc, 

2 rpm, fluid 

502 format (lx, 'Inlet Pressure = ',0pfl0.4,4x, ' kPa' / 

2 lx, 'Exit Pressure = ', 0pfl0.4,4x, ' kPa' / 

3 lx, 'Inlet Radius = ',lpel4.4, ' m' / 

4 lx, 'Exit Radius = *,lpel4.4, 1 m' / 

5 lx, ’Coning Slope = ',lpel4.2, ' m/m' / 

6 lx, 'Conductivity = ',0pfl0.4,4x, ' W/m-K' / 

7 lx, 'Seal Speed = ',0pf8.2, 6x, ' rpm' / 

8 lx, 'Fluid = ',3x,a8 /) 

end if 

pinlet = 1.0d3*pinlet 
pback = 1 . 0d3*pback 

c c 

c Crunch the curves out ... c 

c c 

do 30, i = 1, ncurve 

if ( (choose. eq. 'S' ) .or. (choose. eq. 's') ) then 
omega =3 . 141592653589793d0*rpms (i) /30 . OdO 
write (12, ' (//" RPM = ",fll.4) ') rpms(i) 
write (12, '( " Omega = ",fll.4)') omega 
write (*,'('' RPM(",i2, ") =",fl2.4)') i, rpms(i) 

call match (tinf, rpms(i)) 

if ((both .eq. 'y') .or. (both .eq. 'Y')) then 
call face (tinf, rpms(i)) 
end if 

GlS6 

~ omega =3.141592653589793d0*rpm/30.0d0 
write (12, ' (//" Tinf = ",fll.4, " K"/)') tinfs(i) 
write (*, ' (" Tinf ( " , 12, ' ' ) =",fl2.4)') i, tinfs(i) 
call match (tinfs (i) , rpm) 
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30 


if ((both .eq. ’y’) -or. (both .eq. 'Y*)) then 
call face (tinf s (i) , rpm) 
end if 
end if 
continue 


close (12) 


stop 

end 


c+ 

c 


+c 

c 


subroutine match ( tresvr, rpmsl ) 


The routine which calculates the seal tenperaturea, leakage rates 
forces for various fils, thickness grva. 
operating conditions using the subroutines ALLLIQ and PHAS2L. 


c 
c 
c 
c 
c 

C 4 


c 
c 
c 
c 
c 

he 


cc 


cc 


cc 


illicit logical ( a - z ) 
implicit urvdfefined ( a - z ) 

real *8 r inlet, routlt, cone, conduc 

real* 8 pinlet, pback, rpm, omega, rolinf, tint 

real*8 hinlet, pboil, rtx)il 

real *8 hinupt, hinlot, hinupl, tonlol 

real*8 tseal, visgas, visliq, rolliq, rgas 


integer npntst, nprrtsl 

common /geomty/ r inlet, 
common /operat/ pinlet, 
corrrnon /iterat/ hinlet, 
common /curves/ npntst, 


rout It, cone, conduc 

pback, rpm, omega, rolinf, tinf 

pboil, rtooil # . , . . 

hinupt, hinlot, npntsl, hinupl, hmlol 


corrmon /fluid / tseal, visgas, visliq, rolliq, rgas 


real*8 tresvr, rpmsl, pi, houtlt, rtmax, delh, sat 
integer i 

logical boil, boiled 
external sat 


parameter (pi = 3.141592653589793d0) 


rpm = rpmsl 
tinf = tresvr 
omega = pi*rpm/30.0d0 
boiled = .false, 
rbmax = rout It 

if (tinf .ge. sat(l, pinlet)) then 
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write (\ • (/ " Seal is all-vapor ..." /) ’ ) 
end if 

delh = (hinlol/hinupl) ** (1 .OdO/dble (npntsl - D) 
hinlet = hinupl 


write (12,501) 

501 format (/, ' hrin 

2 ' leakage 

do 10, i = 1, npntsl 

rboil = rout It 


xboil pboil tseal ', 

load'/) 


if (cone .ne. O.OdO) then 

houtlt = hinlet + (routlt - rinlet) *cone 
if (houtlt .le. O.OdO) then 

write (*,' (/ " WARNING: Coning angle specified is , 

2 ' ' too large. ") ') 

1 *^ * ( * * Calculations for film , 

2 1 1 thicknesses less than or equal to " ) 1 ) 

write (*, ' (10x, lpel0.3, ' ' m are skipped ..." /)') 

2 hinlet 

return 
end if 
end if 

if (tinf .ge. sat(l, pinlet) ) then 
ccc 

ccc write (*,'(" All vapor seal ... ")’) 

ccc 

call allvap 


ccc 

ccc 

ccc 


ccc 

ccc 

ccc 


else 

call heat (boil) 

if (boil .or. boiled) then 

write (*, * ( " Two-phase seal ... * * ) ) 

call phas21 (rfcmax, boiled) 
else 

write (*,'(" All liquid seal ... " ) ' ) 


call allliq 
end if 


end if 
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hinlet = hinlet * delh 

10 continue 

return 

end 


c+ 

c 


c 

c 

c 

c 

c+ 


subroutine allliq 

This is an all liquid seal, 
seal speed is exceeded. 


•+c 

c 


But we need to check if the critical 


c 
c 
c 
c 

— +c 


inplicit logical ( a - z ) 
cc implicit undefined ( a - z ) 

real*8 rinlet, routlt, cone, conduc 

real *8 pinlet, pback, rpm, omega, rolinf, tinf 

real *8 hinlet, pboil, rboil 

real *8 tseal, visgas, visliq, rolliq, rgas 

common /geanty/ rinlet, routlt, cone, conduc 

common /operat/ pinlet, pback, rpm, omega, rolinf, tinf 

common /iterat/ hinlet, pboil, rboil 

common /fluid / tseal, visgas, visliq, rolliq, rgas 

real*8 crt, rbold, pi, rerr, load, leak, xboil, wliq, lkliq 
logical boil 
external wliq, lkliq 

parameter (pi = 3.141592653589793d0, rerr = 1.0d-8) 

cr t = pinlet - pback - 0.15d0*rolliq*omega*omega 
2 * (rinlet *rinlet - rout It* rout It) 



c If CRT is negative, the seal is operating at super-critical c 

c speed. 

c 

c 


ccc 
ccc 
ccc 2 
ccc 
10 


2 


if (crt .It. O.OdO) then 

write (*, 1 ( ' ' Crit""ed, Hinlet = ",lpel2.3, 

RPM = ' ',0pfl0.2) ') hinlet, rpm 

rbold = rboil 

rboil = dsqrt (rinlet *rinlet - 20 . OdO* (pinlet - pback) 

/(3.0d0*rolliq*omega*omega) ) 
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call heat (boil) 

if (dabs (l.OdO - rbold/rboil) .gt. rerr) go to 10 

pboil = pback 

load = dabs (wliq(rboil) ) 

+ pback*pi*dabs (rboil*rboil - routlt*routlt) 

xboil = (routlt - rboil) / (routlt - rinlet) 
write (12, 501) hinlet, xboil, tseal, load 
format (2 (lx, lpel2 . 3) , 14x, Opf 10 . 3, 14x, lpel2 .3) 


else 


ccc 

ccc 

ccc 2 
ccc 


502 


write (*,'(" AOKay, Hinlet = ",lpel2.3, 

", RPM = ' ', Opf 10. 2) ') hinlet, rpm 

pboil = pback. 

load = dabs (wliq (routlt)) 

leak = lkliq (routlt) 

write (12, 502) hinlet, tseal, leak, load 
format (lx, lpel2 . 3, 27x, Opf 10 . 3, 2 (lx, lpel2 . 3) ) 


end if 


rboil = routlt 


return 

end 


+c 


subroutine phas21 ( rbmax, boiled ) 

c This seal is likely to be a two-phase seal. But we need to check c 
c if the critcal speed is exceeded. c 

c 

c+ 

implicit logical ( a - z ) 
cc implicit undefined ( a - z ) 

real* 8 rinlet, routlt, cone, conduc 

real*8 pinlet, pback, rpm, omega, rolinf, tinf 

real *8 hinlet, pboil, rboil 

real*8 tseal, visgas, visliq, rolliq, rgas 

common /geomty/ rinlet, routlt, cone, conciuc 

common /operat/ pinlet, pback, rpm, omega, rolinf, tinf 

common /iterat/ hinlet, pboil, rboil 

contnon /fluid / tseal, visgas, visliq, rolliq, rgas 

real*8 rbmax, rerr, lkerr, crt, rbold, load, leak, xboil, 
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2 

3 


dir, taigas, rleft, rright, glkage, check, llkage, sat, 
lkgas, Ddiq, wgas, wliq 
logical boiled, boil 
external lkgas, 13d iq, wgas, wliq 

parameter (rerr = 1 . Od-8, lkerr = 1.0d _ 4) 

crt = pinlet - pback 

2 - 0.15d0*rolliq*cmega*cmega 

3 * (rinlet*rinlet - rboil*rboil) 

ccc 

ccc write (*,'(" CRT = ”,lpel2.3)') crt 
ccc 

if (crt .It. O.OdO) then 

10 rbold = rboil 

rboil = dsqrt (rinlet* rinlet - 20. OdO* (pinlet - pback) 

2 / (3 . OdO*rolliq*omega*omega) ) 

call heat (boil) 

if (dabs (1. OdO - rbold/rboil) .gt. rerr) go to 10 

if (.not .boil) then 
boiled = .false, 
pboil = pback 
load * dabs (wliq (rboil) ) 
xboil = (rout It - rboil) / (routlt - rinlet) 
write (12, 501) hinlet, xboil, tseal, load 
501 format (2 (lx, lpel2 . 3) , 14x, Opf 10 .3, 14x, lpel2 . 3) 

rboil = routlt 
return 
end if 

end if 

if (rinlet .It. routlt) then 
dir = l.OdO 
else 

dir = -l.OdO 
end if 

boiled = .true, 
rleft = rboil 
rright = rinlet 

20 rboil = 0.5d0* (rleft + rright) 

call heat (boil) 
taigas = sat(l, pinlet) 
if (tseal .ge. taigas) then 
rleft = rboil 
go to 20 
end if 

glkage = lkgas () 
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llkage = lkliqO 

check = dir* (glkage - llkage) 

if ( ( (dabs (check/glkage) .gt. lkerr) .or. 

(dabs (check/llkage) .gt. lkerr)) .and. 
(dabsd.OdO - rright/rleft) .gt. rerr)) then 
if (check .gt. O.OdO) then 
rleft = rboil 
else 

rright = rboil 
end if 
go to 20 
end if 


ccc 

ccc 

ccc 2 

ccc 

ccc 

ccc 

ccc 2 

ccc 

ccc 

ccc 

ccc 

ccc 

ccc 

ccc 2 

ccc 

ccc 

ccc 

ccc 


502 


write (*,'(" XBoil = ’\lpel2.3)') 

(routlt - rboil) / (rout It - rinlet) 
write (*, ' (" PBoil = ",lpel2.3)') pboil 


if 


( (dabs (check/glkage) .gt. lkerr) .or. 

(dabs (check/llkage) .gt. lkerr)) then 
if (dabs (check/glkage) .gt. dabs (check/llkage) ) 
check = dabs (check/glkage) 


then 


else 

check = dabs (check/llkage) 


end if 

write (*, ' (/ " 

i r 

write (*, * (llx, 
write (*, ' (llx, 
end if 


WARNING: Leakage rates did not converge", 

to within the specified limit' ') ') 

"at HRIN =' 1 , lpel0.3, ' ' m, ")') hinlet 
" ERR =",fl0.7, /)') check 


load = dabs (wliq( rboil) + wgas (rboil)) 
leak = 0.5d0* (llkage + glkage) 
xboil = (routlt - rboil) / (routlt - rinlet) 
write (12, 502) hinlet, xboil, pboi 1 ' tseal, leak, load 
format (3 (lx, lpel2 . 3) , lx, Opf 10 .3,2 (lx, lpel2 . 3) ) 


return 

end 


c+ 

c 

subroutine allvap 
c 

c This is an all vapor seal, 
c is assumed negligible, 
c 

c+ 




c 

c 

Heat generation by viscous disspation c 

c 




implicit logical ( a - z ) 
implicit undefined ( a — z ) 

real *8 rinlet, routlt, cone, conduc 
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real *8 pinlet, pback, rpm, omega, rolinf, tinf 

real*8 hinlet, pboil, rboil 

real*8 tseal, visgas, visliq, rolliq, rgas 

cortnon /geanty/ rinlet, routlt, cone, conduc 

common /operat/ pinlet, pback, rpm, omega, rolinf, tinf 

common /iterat/ hinlet, pboil, rboil 

common /fluid / tseal, visgas, visliq, rolliq, rgas 

real *8 load, leak, sat, wgas, lkgas 
external sat, wgas, lkgas 

tseal = tinf 
visgas = sat (8, tseal) 
rboil = rinlet 
pboil = pinlet 
load = wgas (rinlet) 
leak = lkgas (rinlet) 

write (12, 501) hinlet, tseal, leak, load 
501 format (lx, lpel2.3,27x, OpflO.3,2 (lx, lpel2.3) ) 

return 

end 


c+ 

c c 

subroutine heat ( boil ) 

c c 

c Subroutine that calculates, for a given boiling interface location, c 
c the temperature at the boiling interface using the semi-infinite c 
c solid heat transfer model. The leakage flow is assumed to be c 

c quasi-isothermal at this terrperature . The pressure at the boiling c 

c interface and the saturation fluid properties are then obtained c 

c from the steam table ( real*8 function SAT ) . c 

c c 

c+ +c 


implicit logical ( a - z ) 
cc implicit undefined ( a - z ) 

real *8 rinlet, routlt, cone, conduc 

real*8 pinlet, pback, rpm, omega, rolinf, tinf 

real*8 hinlet, pboil, rboil 

real*8 tseal, visgas, visliq, rolliq, rgas 

real*8 gamma, cpgas, tcrit, vcrit, pcrit, psmall 

common /geanty/ rinlet, routlt, cone, conduc 

common /operat/ pinlet, pback, rpm, omega, rolinf, tinf 

common /iterat/ hinlet, pboil, rboil 

common /fluid / tseal, visgas, visliq, rolliq, rgas 

common /gas / gamma, cpgas, tcrit, vcrit, pcrit, psmall 
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real*8 ttemp, intgrl, deltt, terr, t small, pi, blkpt 1 , blkpt 2 , 
blkpt3, taigas, ttrial, termht, sat, gss4, hcondc 
real*8 intgll, intgl2, intgl3, intgl4 
logical boil 

external sat, gss4, hcondc 

parameter (pi=3.141592653589793d0, terr=1.0d-2, tsmall=10.0d0) 


c Since the integrand function has a singularity at R - RBOIL, an c 
c open interval integration scheme (4 points Gauss-Legendre c 
c integration over 4 sub-intervals) is chosen. c 


blkpt 1 = rboil + 0 . ldO* (rinlet - rboil) 
blkpt2 = rboil + O.OldO* (rinlet - rboil) 
blkpt 3 = rboil + O.OOldO* (rinlet - rboil) 
intgrl = gss4 (rinlet, blkptl, hcondc) 

2 + gss4 (blkptl, blkpt2, hcondc) 

3 + gss4 (blkpt 2, blkpt3, hcondc) 

4 + gss4 (blkpt3, rboil, hcondc) 


cc 

cc 

cc 

cc 

cc 

cc 

cc 

cc 

cc 

cc 

cc 

cc 

cc 

cc 


cc 


cc 

cc 


cc 


intgll = gss4 (rinlet, blkptl, hcondc) 

write (*,’(' ’INT(LCW =",lpelO. 3, ' ', UP = ,lpel0.3, ) - , 

2 lpel0.3) 1 ) rinlet, blkptl, intgll 

intgl2 = gss4 (blkptl, blkpt2, hcondc) 

write (*,'(" INT (LCW =' ' , lpel0.3, 1 ' , UP = ,lpel0.3, ) - , 

2 lpel0.3) ') blkptl, blkpt2, intgl2 

intgl3 = gss4(blkpt2, blkpt3, hcondc) 

write (*,'('' INT (LOW =' ' , lpel0.3, 1 1 , UP = ,lpel0.3, ) - , 

2 lpel0.3) ') blkpt2, blkpt3, intgl3 

intgl4 = gss4 (blkpt3, rboil, hcondc) ,, 

write (*, * (' ' INT(LCW =' ',lpel0.3, ' ', UP = ',lpel0.3, ') = , 

2 lpel0.3) ') blkpt 3, rboil, intgl4 

write (*, ' (lx) ’) 

intgrl = intgll + intgl2 + intgl3 + intgl4 


if (rboil .gt. rinlet) then 
intgrl = -l.OdO* intgrl 
end if 

write (*, ' (/" INTGRL = ",f20.8)') intgrl 


2 


taigas = sat(l, pinlet) 
visliq = sat (7, taigas) 

write (*,'(" VisLiq (TAlGas =",f9.3, ' ) = 
taigas, visliq 

termht = omega* omega* intgrl/ (pi*conduc) 
deltt = visliq*termht 
ttenp = tinf + deltt 


, f20.8) ') 


if (tterp .gt. taigas) then 

write (*, ' ( ' ' Seal will be all gas — * * ) 1 ) 
tseal = taigas 
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boil = .true, 
return 
end if 


20 


cc 

cc 

cc 

cc 


visliq — sat (7 , ttenp) 
deltt = visliq*termht 
ttrial = tinf + deltt 

if (dabs(ttenp - ttrial) .gt. terr) then 
if (ttrial .gt. taigas) then 

ttenp = 0 . 5d0* (taigas + ttemp) 
else 

ttenp = ttrial 
end if 
go to 20 
end if 


tseal = ttrial 
pboil = sat (6, tseal) 
visgas = sat (8, tseal) 
rolliq = l.OdO/sat (10, tseal) 
if (pboil .gt. pback) then 
boil = .true, 
else 

boil = .false, 
end if 

write (*, ' (" Tseal = ,, ,f20.8)') 
write (*,'(" Pboil = '\f20.8)') 
write (*,'('' VisGas = ,, ,f20.8) ) 
write (*,'(' 1 RolLiq = ",f20.8)') 


tseal 

pboil/1000 . OdO 

visgas 

rolliq 


return 

end 


c+ 

Q 

real *8 function hcondc ( radius ) 

c The integrand function of the heat transfer model, 
c 

c+ 


cc 


cc 

cc 


cc 


implicit logical ( a - z ) 
implicit undefined ( a - z ) 

real *8 rinlet, routlt, cone, conduc 

real*8 pinlet, pback, rpm, omega, rolinf, tinf 

real* 8 hinlet, pboil, rboil 

real*8 tseal, visgas, visliq, rolliq, rgas 


ommon /geoirty/ rinlet, routlt, cone, conduc 

omnon /operat/ pinlet, pback, rpm, omega, rolinf, tinf 

omnon /iterat/ hinlet, pboil, rboil 
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cctrmon /fluid / tseal, visgas, vislig, rolliq, rgas 


real *8 radius, ellipl 
external ellipl 


cc 

cc 


if (radius .It. rboil) then 

hcondc = radius*radius*radius*ellipl (radius/rboil) 

/ (rboil* (hinlet + (radius - rinlet) *cone) ) 

6lSG 

hcondc = radius*radius*ellipl (rboil/ radius) 

/(hinlet + (radius - rinlet) *cone) 

end if 

write (*, ' (" HCONDC (r = "iflO.B, ") =",f20.8)') radius, 
hcondc 


return 

end 


c+- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c+- 


real* 8 function ellipl ( x ) 

A double precision function that evaluates the complete elliptical 
integral of the 1st kind using the series approximation given in 
the "Handbook of Mathematical Functions." 

The absolute error of this approximation is less than 3.0e-5. 

Reference : Abramowitz, M. , and Stegun, I. A., 

"Handbook of Mathematical Functions", pp 591, 

Equation 17.3.34, Ninth Printing, Dover Publications 

Note that the relationship between the parameter Ml as defined 
in equation 17.3.34 and the input parameter of this function, X, 
is i 

Ml = 1.0 - X**2 


-+c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

-+c 


implicit logical ( a - z ) 
cc implicit undefined ( a - z ) 

real *8 x, ml, m2, aO, al, a2, bO, bl, b2 

parameter ( aO = 1.38629 44 dO, bO = 0.50000 00 dO, 

2 a l = 0.11197 23 dO, bl = 0.12134 78 dO, 

3 a2 = 0.07252 96 dO, b2 = 0.02887 29 dO) 


ml = l.OdO - x * x 
m2 = ml * ml 

ellipl = (aO + al * ml + a2 * m2) 

- (bO + bl * ml + b2 * m2)*dlog(ml) 
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return 

end 


C + “ 

c 

subroutine rsatdt ( fluid ) 

c Subroutine that reads in the spline coefficients for the steam 
c tables from the saturation fluid properties data files, 
c 

c+ 



inplicit logical ( a - z ) 
implicit undefined ( a - z ) 

real *8 gamma, cpgas, tcrit, vcrit, pent, psmall 
real *8 satdat, ddatdp, visdat, pvisc, tvisc, tmin, tmax, 
2 prnin, penax 

real *8 tseal, visgas, visliq, rolliq, rgas 
integer ndata, jprops, irows, jcols, i, j ,k 
character* 8 fluid 


2 

3 


common /gas / gamma, cpgas, tent, vent, pent, psnall 
■'ormvon /thermo/ satdat (11,50,4) , ndata, jprops, ddatdp (5, 60, 3) , 
visdat (65, 15) , pvisc (15), tvisc (65), irows, 
jcols, tmin, tmax, pnin, pmax 
oommon /fluid / tseal, visgas, visliq, rolliq, rgas 


jprops = 1 

if (fluid .eq. 'water ') then 

open (10, f i le= ' water . dat ' , status='old') 
rgas — 461.51d0 

else if (fluid .eq. 'hydrogen') then 

open (10, f ile=' hydrogen. dat ' , status— old ) 
rgas = 4124.289d0 

else if (fluid .eq. 'oxygen ') then 

open (10, f ile= ' oxygen . dat ' , status='old’) 
rgas = 259.832d0 

else if (fluid .eq. 'nitrogen') then 

open (10, file= 'nitrogen.dat ', status-'old' ) 
rgas = 296.798d0 
0lSG 

write (*, ' ( ' ' Program cannot handle this fluid ... ' ' ) ) 
stop 
end if 

read(10,701) gamma, cpgas, tcrit, vcrit, perit, psmall 

read (10, 702) ndata .. 

read (10,703) ( ( (satdat (i, j, k) , i=l, 11) , 3 = 1 , ndata) , k-1, 4 
read (10, 703) ( ( (ddatdp (i, j,k) , i=l, 5) , j=l, ndata) ,k-l,3) 
read (10, 702) irows, jcols 
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10 


read (10, 701) tmin, tmax, pmin, pmax 
do 10, j = 1, jcols 
read(10,701) pvisc(j) 
do 20, i = 1, irows 
20 read (10, 701) tvisc(i) 

do 30, j = 1, jcols 
do 30, i = 1, irows 
30 read(10,*) visdat(i,j) 

close (10) 

return 

701 format (bn, el3. 6) 

702 format (bn, il3) 

703 format (bn, 6 (lpel3. 6) ) 

end 


c+- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c+ 


■+c 

c 


real*8 function sat(i, x) 

This function calculates saturation properties as functions of 
either pressure, temperature, or saturated liquid enthalpy. The 
data is stored in the array SATDAT (I, J,K) where I denotes the 
the fluid property as follows: 


1 

2 

3 

4 


Tsat 

Vf 

Vfg 

If 


(Psat) 

(Psat) 

(Psat) 

(Psat) 


5 

6 

7 

8 


Ifg 

Psat 

VISCf 

VISCg 


(Psat) 

(Tsat) 

(Tsat) 

(Tsat) 


9 

10 

11 


If 

Vf 

Tsat 


(Tsat) 

(Tsat) 

(If) 


The index K takes on values 1 through 4 where K - 1 denotes the 


J. K = 2,3,4 


propery for the particular 
spline coefficient for the curve of property type I 
between point J and J + 1. 


contains the cubic 
(e.g. If) 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

-+c 


implicit real*8 (a-h,o-z) 

common /gas / gamma, cpgas,tcrit,vcrit,pcrit,psmall 
common /thermo/ satdat (11, 50, 4) ,ndata, jprops,ddatcip(5, 60, 3) , 

2 visdat (65, 15) , pvisc (15) , tvisc (65) , irows, jcols, 

3 tmin, tmax, pmin, pmax 


dimension ipropx(ll) 

data ipropx /6, 6, 6, 6, 6, 1,1, 1, 1, 1,4/ 


c 

c Remember what property 
c I. 

c 




is used as the independent variable for this c 

c 
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o o o 


ix = iprqpx(i) 




^ in j 31x1 c 

c J + 1. 

if ((x .ge. satdat (ix, jprops ,1)) .and. 

2 (x .le. satdat (ix, jprcps+1, D ) ) go to 20 

Find the correct J for this property^ c 

jmove = int (dsign(1.0d0, x - satdat (ix, jprops, 1) ) ) 

100 jprops = jprops + jmove 

if ( (x .It. satdat ( ix, jprcps ,1)) -or. 

2 (x .gt. satdat (ix, jprcps+1, 1))) go to 10 

c Calculate the property from toe^cubic^splitie^ 

” , 


return 

end 


c+ 

c 


real* 8 function wliq ( rhoil ) 




C 

c 
c 
c 

c+ 


+c 

c 

c 

c 

c 

c 

+c 


implicit logical ( a - z ) 
implicit undefined ( a - z ) 

real *8 rboil, rortibrg, tprliq 
real*8 rinlet, routlt, cone, conduc 
external rombrg, tprliq 

cannon /geomty/ rinlet, routlt, cone, conduc 
w liq = ranbrg (rinlet, rboil, tprliq) 
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return 

end 


c+- 

c 

c 

c 

c 

c 

c 

c 

c 

c+- 


real*8 function tprliq ( radius ) 

Function that calculates the value of : 2.0 * pi * r * P (r) 

for the all-liquid portion of the seal. . ^ 

To find the seal opening force frcm the all liquid ° 

seal, this function is integrated with respect to RADIUS frcm 
the inlet to the boiling interface location. 


■+c 

c 

c 

c 

c 

c 

c 

c 

c 

-+c 


cc 


iirplicit logical ( a - z ) 
implicit undefined ( a - z ) 


real* 8 radius, press, twopi, pliq 
integer ncalls 
external pliq 

parameter (twopi = 6.283l85307179586d0) 
data ncalls /0/ 


cc 

cc501 
cc 2 


ncalls = ncalls + 1 
press = pliq (radius) 
tprliq = twopi * radius * press 
write (*,501) ncalls, radius, press, tprliq 
formatC Call # f ,i3, P (r = * ,fl0.8, f ) - 
23x, ’2*pi*P = ' ,f20.6) 


,f20.6 


/ 


return 

end 


c+ 

c v 

real *8 function pliq ( radius ) 

c Function that calculates the pressure at RADIUS for the all 
c liquid portion of the seal. 

c 



implicit logical ( a - z ) 
implicit undefined ( a - z ) 

real *8 r inlet, rout It, cone, conduc 

real *8 pinlet, pback, rpm, omega, rolinf, tinf 

real *8 hinlet, pboil, rboil 

real*8 tseal, visgas, visliq, rolliq, rgas 
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2 

3 

4 


2 

3 

2 

3 

4 

2 

3 

4 


corrmon /geomty/ r inlet, rout It, cone, conduc 

corrrron /operat/ pinlet, pback, rpm, omega, rolinf, tint 

corrmon /iterat/ hinlet, pboil, rboil 

common /fluid / tseal, visgas, visliq, rolliq, rgas 

real^8 radius, terml, centrf, h, ho, hboil, numerl, nuner2, 
denom 

if (cone .eq. O.OdO) then 
terml = (pinlet - pboil 

- 0.15d0*rolliq*cmega*cmega 

* (rinlet*rinlet - rboil*rboil) ) 
★diog (radius/rinlet) /dlog (rinlet/rboil) 

centrf = 0 . 15 d 0 *rolliq*cmega*cmega . . 

* (radius*radius - r inlet * r inlet ) 

pliq = pinlet + terml + centrf 

else , , 

h = hinlet + (radius - r inlet) *cone 
hboil = hinlet + (rboil - rinlet)*cone 
ho = hinlet - rinlet*cone 
numerl = pinlet - pboil 

- 0 . 15d0*rolliq*cmega*cmega 

* (rinlet*rinlet - rboil*rboil) 

numer2 = dlog ( (radius*hinlet) / (rinlet*h) ) 

+ ho* (l.OdO/h - 1 . OdO /hinlet) 

+ 0 . 5d0*ho*ho* (1 . OdO/ (h*h) 

- l.OdO/ (hinlet *hinlet ) ) 

denom = dlog ( (rinlet*hboil) / (rboil*^let) ) 

+ ho* (1 . OdO /hinlet - 1. OdO /hboil) 

+ 0 . 5d0*ho*ho* (1 . OdO/ (hinlet *hinlet) 

- l . OdO/ (hboil*hboil) ) 

centrf = 0 . 15 d 0 *rolliq*cmega*croega , 

* (radius*radius - rinlet*rinlet) 

pliq = pinlet + numerl *numer2 /denom + centrf 

end if 


return 

end 




c+ - c 

c 

real* 8 function lkliq ( ) c 

c Function that calculates the mass leakage rate from the liquid seal c 
c equations. c 

C 


c 


implicit logical ( a - z ) 
implicit undefined ( a - z ) 
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real *8 r inlet, rout It, cone, conduc 

real*8 pinlet, pback, rpm, anega, rolinf, tinf 

real *8 hinlet, pboil, rboil 

real *8 tseal, visgas, visliq, rolliq, rgas 

cannon /geomty/ rinlet, routlt, cone, conduc 

common /qperat/ pinlet, pback, rpm, omega, rolinf, tinf 

common /iterat/ hinlet, pboil, rboil 

conmon /fluid / tseal, visgas, visliq, rolliq, rgas 

real *8 pi, ho, hboil, numer, dencm 

parameter (pi = 3.141592653589793d0) 


CC 


write (*,501) pi, rolliq, hinlet, omega. 

rinlet. 

cc 

2 

pboil 



cc501 

format ( ' Pi = ' , f 20 . 8 / 

' RolLiq = 

',f20.8 

cc 

2 

' Hinlet = ' , f20.8 / 

' Omega = 

',f20.8 

cc 

3 

• Rinlet = ' ,f20.8 / 

• Rboil = 

',f20.8 

cc 

4 

• Pinlet = ' ,f20.8 / 

' Pboil = 

',f20.8 


rboil, pinlet, 

/ 

/ 

/ 

) 


cc 

cc 

cc 


cc 

cc 

cc 


cc 

cc 

cc 


2 

3 

4 


2 

3 

4 

2 

3 

4 

5 


if (cone .eq. O.OdO) then 

numer = rolliq*pi*hinlet*: ?u*hinlet 

* (0.15d0*rol^iq*omega*omega 

* (rinlet *rinlet - rboil*rboil) 

- (pinlet - pboil) ) 

dencm = 6.0d0*visliq*dlog(rinlet/rboil) 
write (*, ' (" Numer = ,, ,f20.8) 1 ) numer 
write (*, '(" Dencm = ",f20.8)') dencm 
lkliq = numer/denom 

write (*,'(" LkLiq = ,, ,f20.8)') lkliq 
else 

ho = hinlet - rinlet*cone 

hboil = ho + rboil*cone 

write (*,'(" Hinlet = ",f20.8)') hinlet 

write (*, '(" Ho = ",f20.8)') ho 

write (*, •(" Hboil = ",f20.8)') hboil 

numer = rol 1 iq*pi *ho*ho*ho 

* (0.15d0*rolliq*omega*omega 

* (rinlet*rinlet - rboil*fboil) 

- (pinlet - pboil) ) 


dencm = 6.0d0*visliq 

* (dlog (rinlet /rboil) - dlog (hinlet/hboil) 

+ ho* (l.OdO /hinlet - l.OdO/hboil) 

+ 0 . 5d0*ho*ho* ( 1 . OdO/ (hinlet *hinlet) 

- 1 . OdO/ (hboil*hboil) ) ) 
write (*, ' (" Numer = ,, ,f20.8)') numer 
write (*, ' (" Denom = ' , ,f20.8)') denom 
lkliq = numer /dencm 

write (*,'(" LkLiq = ",f20.8)') lkliq 
end if 


return 

end 
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c+ 

c 


real *8 function wgas ( rboil ) 


c 
c 
c 
c 

c 


Function that calculates the opening force from the all vapor 
portion of the seal given the boiling interface location. 


— +c 
c 

c 

c 

c 

c 


implicit logical ( a - z ) 
implicit undefined ( a - z ) 

real *8 rboil, roribrg, tprgas 
real*8 rinlet, routlt, cone, conduc 
external rombrg, tprgas 

common /geomty/ r inlet, rout It, cone, conduc 

wgas = rcrribrg (rboil, rout It, tprgas) 

return 

end 


c+ 

c 


+c 

c 


real *8 function tprgas ( radius ) 


c 
c 
c 
c 
c 
c 
c 

c+ 


2.0 * pi 


P (r) 


Function that calculates the value of 

To°f ind. e the 1 sS 1° open ing°f orce from the 0 P °Sj?Ss° f f^ 

seal, this function is integrated with respect to RADIUS 
the inlet to the boiling interface location. 


c 
c 
c 
c 
c 
c 
c 




cc 


implicit logical ( a - z ) 
irrplicit undefined ( a - z ) 


real* 8 radius, twopi, press, pgas 
integer ncalls 
external pgas 

parameter (twopi = 6.283185307179586d0) 
data ncalls /0/ 


ncalls = ncalls + 1 

press — pgas (radius) 

tprgas = twopi * radius * press 

write (*,501) ncalls, radius, press, tprgas 
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cc501 
cc 2 


formate Call #',i3, ' , P(r - ’ ,fl0.8, ') -',f20.6 / 
23x, '2*pi*P = ’,f20.6) 


return 

end 


c+ 

c 

real*8 function pgas ( radius ) 
c 

c Function that calculates the pressure at RADIUS 
c vapor portion of the seal, 
c 

c+ 


+c 

c 

c 

for the all- c 

c 

c 

+c 


implicit logical ( a - z ) 
cc implicit undefined ( a - z ) 


cc 


real *8 rinlet, rout It, cone, conduc 

real *8 pinlet, pback, rpm, omega, rolinf, tinf 

real*8 hinlet, pboil, rboil 

real*8 tseal, visgas, visliq, rolliq, rgas 


common /geomty/ rinlet, rout It, cone, conduc 

common /operat/ pinlet, pback, rpm, omega, rolinf, tinf 

common /iterat / hinlet, pboil, rboil 

common /fluid / tseal, visgas, visliq, rolliq, rgas 

real* 8 radius, terml, ho, h, hboil, houtlt, numerl, numer2, 
2 denom 


if (cone .eq. O.OdO) then 

terml = (pboil*pboil - pback*pback) *dlog(radius/routlt) 

2 /dlog(rboil/routlt) 

pgas = dsqrt (pback *pback + terml) 
else 

ho = hinlet - rinlet* cone 

h = ho + radius* cone 

hboil = ho + rboil*cone 

houtlt = ho + routlt*cone 

numerl = pboil *pboil - pback *pback 

numer2 = dlog (radius*houtlt/ (routlt*h) ) 

2 + ho*(1.0d0/h - l.OdO/houtlt) 

3 + 0 . 5d0*ho*ho* (1 . OdO/ (h*h) 

4 - 1 . OdO/ (houtlt*houtlt) ) 
denom = dlog (rboil*houtlt/ (routlt*hboil) ) 

2 + ho* (1. OdO /hboil - 1. OdO /houtlt) 

3 + 0 . 5d0*ho*ho* (1 . OdO/ (hboil*hboil) 

4 - 1. OdO/ (houtlt *hout It) ) 
pgas = dsqrt (pback *pback + numerl *numer2/denan) 

end if 
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return 

end 


c+ 

c 

real *8 function Ikgas ( ) 

c Function that calculates the mass leakage rate from the vapor 
c equations. 

c 

c+ 



c 

seal c 
c 
c 

+c 


cc 


implicit logical ( a - z ) 
irrplicit undefined ( a - z ) 

real *8 r inlet, routlt, cone, conduc _ 

real *8 pinlet, pback, r F n *' i ane 9 a ' rolinf, tint 

real* 8 hinlet, pboil, rboil 

real *8 tseal, visgas, visliq, rolliq, rgas 


common /geomty/ rinlet, routlt, cone, conduc . _ 

common /operat/ pinlet, pback, rpm, omega, rolinf, tint 
common /iterat/ hinlet, pboil, rboil 
conmon /fluid / tseal, visgas, visliq, rolliq, rgas 


real *8 pi, ho, hboil, houtlt, numer, denom 


parameter (pi=3.141592653589793d0) 


cc 

cc 2 
cc501 
cc 2 


cc 

cc 

cc 


3 

4 

5 


write (*,501) pi, pback, pboil, hinlet, cone, visgas, 
tseal, rboil, routlt 


format ( ' Pi 

' Pboil 
' Cone 
' Rgas 
' Pboil 


1 , f 20 . 8 / ' Pback = ',f20.8 / 

1 ,f20.8 / ' Hinlet = ',f20.8 / 

' , f20.8 / ' VisGas = ',f20.8 / 

' ,f20.8 / ' Tseal = \f20.8 / 

' ,f20.8 / ' Routlt = 1 ,f20.8 ) 


rgas, 


cc 

cc 

cc 


cc 

cc 

cc 


2 


^ numer = pi* (pback*pback - pboil*pboil) *hinlet*hinlet 

*hinlet 


dencxn = I 2 . 0 d 0 *visgas*rgas*tseal*dlog(rboil/routlt) 

write (*, ' (" Numer = 1 \f20.8)’) numer 
write <V<" Dencm * "#*20.8)') denom 

Ikgas = numer/denom 

write (*,'('' LkGas = ",f20.8)'> Ikgas 
else 


ho = hinlet - rinlet*cone 
hboil = ho + rboil* cone 
houtlt = ho + rout It* cone 
write (*, 1 ( f 1 Ho = ,f ,f2Q.8)’) 
write (*, 1 ( 1 1 Hboil — * 1 / f20 .8) ) 
write (*, * (' * Houtlt — ,, r f20.8)') 


ho 

hboil 

houtlt 
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numer = pi* (pfoack*pback - pboil*pboil) *ho*ho*ho 
ctenom = 12 . OdO*visgas*rgas*tseal 

* (dlog (riboil/routlt) - dlog(hboil/houtlt) 

+ ho* (l.OdO/hboil - l.OdO/houtlt) 

+ 0 . 5d0*ho*ho* (l.OdO/ (hboil*hboil) 

- 1 . OdO/ (houtlt*houtlt) ) ) 
write (*, ' (" Numer = ' ',f20.8)') numer 
write (*,'(" Dencm = ",f20.8)') denom 
lkgas = numer /dencm 

write (*, ' (" LkGas = ”,f20.8)’) lkgas 
end if 


return 

end 


c+ 

c 

real* 8 function gss4 ( up, low, func ) 

0 

c This routine calculates the integral of FUNC fran I£W to 
c using 4 points Gauss “Legendre Quadratures • 
c 

c+ 


+c 

c 

c 

UP c 

c 
c 

+c 


irrplicit logical ( a - z ) 
cc inplicit undefined ( a - z ) 

real *8 up, low, xmid, hdelx, xpt, gl, g2, g3, g4, sum, func 
real*8 xi, wt 
dimension xi(2), wt(2) 
external func 

data xi /0. 33998 10435 84856 dO, 0.86113 63115 94053 d0/ 

data wt /0. 65214 51548 62546 dO, 0.34785 48451 37454 dO/ 

xmid = 0.5d0*(up + low) 
hdelx = 0.5d0*(up - low) 

xpt = xmid + hdelx * xi (1) 

gl = func (xpt) 

xpt = xmid - hdelx * xi(l) 

g2 = func (xpt) 

xpt = xmid + hdelx * xi(2) 

g3 = func (xpt) 

xpt = xmid - hdelx * xi (2) 

g4 = func (xpt) 

sum = wt(l) * (gl + g2) + wt(2) * (g3 + g4) 
gss4 = hdelx * sum 


return 
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end 


c+ 

c 


+c 

c 


real* 8 function rcmbrg ( up, low, func ) 

Ro mber g Integration routine. Modified from "Numerical Recipes." 

Reference : Press, w. H., Flannery, B. P., Teukolsky, S. A., and 

Vetter ling, W. T., # ^ 

Numerical Recipes, The Art of Scientific Confuting 
( 1st ed., Cambridge University Press, 1986 ) 
V pp 114-115, 


c 
c 
c 
c 
c 
c 
c 
c 
c 

c+ 


c 
c 
c 
c 
c 
c 
c 
c 
c 

+c 


irrplicit logical ( a - z ) 
implicit undefined ( a - z ) 

real *8 up, low, errmax, s, h, ss, dss, func 
integer nmax, nmaxpl, k, kml, j 
external func 

parameter (errmax=1.0d-8, nmax=20, nmaxpl=21, k=5, kml=4) 
dimension s (nmaxpl), h (nmaxpl) 


10 


h(l) = l.OdO 
do 10, j = 1, nmax 

call trapzd(func, low, up, s(;j), D) 
if (j -ge. k) then 

call polint (h(j-kml) , s(j-kml), k, O.OdO, 
if (dabs (dss) .It. errmax*dabs (ss) ) then 
rombrg = ss 
return 
end if 
end if 

s ( j+1) = s(j) 

h( j+1) = 0.25d0*h(j) 

continue 


ss ,dss) 


write (*,'(' 1 Too many steps ... ' 


stop 

end 


c+ 

c , 

subroutine trapzd (func, low, up, s, n) 

c 



c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c+ 


This routine computes the N’th stage of refinement ofan extended 
trapezoidal rule. Modified from "Numerical Recipes. 

Reference : Press, W. H., Flannery, B. P., Teukolsky, S. A., and 

Numerical Recipes, The Art of Scientific Ccxrputing 
( 1st ed. , Cartbridge University Press, 1986 ) 

pp 111. 


c 

c 

c 

c 

c 

c 

c 

c 

c 

+c 


implicit logical ( a - z ) 

;c implicit undefined ( a - z ) 

real *8 func, low, up, s, tnm, del, x, sum 
integer n, it, j 
external func 

if (n.eq.l) then 

s = 0.5d0*(up - low) * (func (low) + func (up) ) 
it = 1 
else 

tnm = dble (it) 
del = (up - low) /tnm 
x = low + 0.5d0*del 
sum = O.OdO 
do 10, j = 1, it 
sum = sum + func (x) 
x = x + del 

10 continue 

s = 0.5d0* (s + (up - low) *sum/tnm) 

it = 2*it 
end if 


return 

end 


This routine does polynomial interpolation or extrapolations. 
Modified from "Numerical Recipes." 


c+ 

c 

subroutine polint ( xa, ya, n, 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c+ 


+c 

c 


X, y, dy ) 


Reference 


Press, W. H., Flannery, B. P., Teukolsky, S. A., and 

Vetterling, W. T., , _ 

Numerical Recipes, The Art of Scientific Conputing 
1st ed., Cambridge University Press, 1986 ) 

pp 82. 


( 


c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

+c 
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inplicit logical ( a - z ) 
cc inplicit undefined ( a - z ) 

real*8 xa, ya, x, y, dy, dif, dift, c, d, den ,ho, hp, w 
integer nmax, n, ns, i, m 

parameter (nmax=10) 

dimension xa (n) , ya (n) , c (nmax) , d (nmax) 
ns = 1 

dif = dabs (x - xa (1) ) 
do 10, i = 1, n 

dift = dabs(x - xa(i)) 
if (dift. It. dif) then 
ns = i 
dif = dift 
end if 
c(i) = ya(i) 
d(i) = ya(i) 

10 continue 

y = ya(ns) 
ns = ns - 1 
do 30, m = 1, n-1 
do 20, i = 1, n-tn 
ho = xa(i) - x 
Ip = xa(i+m) - x 
w = c(i+l) - d(i) 
den = ho - hp 
if (den .eq. (O.OdO)) then 

write (*,'(" XA" "s are identical ' ' ) ' ) 

stop 
end if 
den = w/den 
d(i) = Ip * den 
c(i) = ho * den 
20 continue 

if ( (2*ns) .It. (n-m) ) then 
dy = c(ns + 1) 
else 

dy = d(ns) 
ns = ns - 1 
end if 
y = y + dy 
30 continue 

return 

end 


subroutine face (tinf, rpn) 
inplicit real *8 (a-h,o-z) 
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carmen 

cannon 

? 

cannon 

? 

cannon 


$ 

$ 


cannon 

cannon 


/gas / 
/geomry/ 

/flow / 

/thermo/ 


/var / 
/stress/ 


gairtna, cpgas, tcrit, verit, pent 

hrrean, radius, rin, rout, r inlet, dhdr, nspace, 

delta, epsi, pi, rmean 

flux, cmega, pinfty, pback, iinfty, txnfty, cho e, 
cf lux, pchoke, closs, signum, 

satdat(ll,50,4), ndata, jprpps, ddatdp(5,60, 3 ) , 
visdat (65,15) , pvisc(15), tvisc(65), irows, Dcols, 
tmin, tmax, pmin, pmax 
p (120) , h(120) 
i-anzr. tauzth# v, mu 


c Contnon blocks to connect the variables free. Lau's laminar program c 
c to Beatty's turbulent program. c 
c 


common /geernty/ srin, srout, scone, sco ™*^ . f 

coramn /operat/ spinlt, spback, srpn, sonega, srol, strnf 
cannon /fluid / stseal, svisg, svisl, sroll, srgas 
cannon /turbin/ nsspac, scloss, sepsi, serror 

/curves/ npntst, hinupt, hinlot, nfantsl, hinupl, hinlol 


comnon 


c 


c+ 



real* 8 mu, iinfty, isat 


c+ - c 

c Put the values of each variable in Lau's program into the corres- c 
c ponding ones in beatty 1 s program . c 


icase = 1 

if (srin. It. srout) then 
rin = srin 
rout = srout 
pin = spinlt 
pout = spback 
else 

rin = srout 
rout = srin 
pin = spback 
pout = spinlt 
end if 
dhdr = scone 
nspace = nsspac 
closs = scloss 
epsi = sepsi 
error = serror 
tinfty = tinf 
iplot = 0 
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c 

c+- 


c 

•+c 


pi = 4 . OdO*datan (1 . OdO) 
omega = pi*rpm/30. OdO 

c Specify finite difference grid. By making each successive step 

c smaller than the last by a set factor, high gradients near the 
c exit may be better resolved. 

sum = (l.OdO - epsi**nspace) / (l.OdO - epsi) 
delta = (rout - rin) /sum 
rmean = (rout + rin) /2. OdO 

c The error bound for the iteration to match the pressure boundary 

c condition (error) has already been set in the input file. Set the 
c error bound for the critical mass velocity and exit pressure, 
bound. 


bound = 1.0d-04 


c+ +c 

c c 

c Modifications made to this subroutine (which was originally a c 

c separate program) so that it calculates the seal performance for c 

c a range of film thicknesses instead of at only one thickness as c 

c in the original program. This is done to make the task of 

c plotting stiffness curves a little easier, 
c c 


hstep = (hinupt/hinlot) ** (1 . OdO/dble (npntst - 1) ) 
hrin = hinupt 


c c 

c+ +c 

c Determine the direction of flow through the seal. 


if (pin. eq. pout) check = O.OdO 

if (pin. ne. pout) check = (pin - pout) /dabs (pin - pout) 
if (check) 12,13,13 

12 pinfty = pout 
pback = pin 
rinlet = rout 
signum = -l.OdO 
go to 14 

13 pinfty = pin 
pback = pout 
rinlet = rin 
signum = l.OdO 
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go to 14 


14 delta = signum*delta 


c+ 

c 

c The big do loop 
c 


■+c 

c 

c 

c 


write (12, 510) 

510 format (/,' hr in Reyc Reyr qexit , 

2 ' leakage load'/) 

do 100 iloops = 1, npntst 

hmean = hrin + 0.5d0*dhdr* (rout - rin) 
iterat = 1 


c FLOOR is hard coded in this subroutine. It MAY NOT be appropriate c 
c for all conditions. If you found that the turbulent part of this c 
c program doesn't run properly, you may have to adjust the c 
c values in the following statements. c 


if (hmean .le. (1.0d-4)) then 
floor = 1.0d-4 
else 

floor = 1.0d-2 
end if 



write (*,'('' Hmean = ",lpel2.3, ' Iloops = ",i4)') 
2 hmean, iloops 



Q 

c Flush all buffered output to save the already computed stuff c 

c into the output file in case this routine enters an infinite loop c 

c and needs to be killed. c 

c Does not work when running under DOS, of course. c 

c 

cc call flush (12) 



Initialize the check for choking. 


choke = -l.OdO 

Calculate the specific enthalpy of the fluid in the reservoir. 
The sub-cooled liquid there is assumed to be incompressible. 
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volume = sat (10,tinfty) 
deltap = pinfty - sat (6,tinfty) 
isat = sat (9,tinfty) 
iinfty = isat + volume*deltap 

Find the lowest back pressure for which a physically realize 
able matching exit pressure may be found. 

fluxl = floor 
flux = fluxl 
call shoot (pexit) 

if (choke. ne. (l.dO) ) go to 21 

write (*,9) 

return 


21 fluxr = fluxl 

22 flux = fluxr 
call shoot (pexit) 
fluxr = 3.162d0*fluxr 

if (choke. ne. (1 .dO) ) go to 22 
choke = -l.OdO 

c Find the mass flux for incipient choking by the bisection 

c method. 

fluxr = fluxr/3.l62d0 

23 flux = fluxl 
cc 

iterat = iterat + 1 
cc 

call shoot (pexit) 
a = choke 
choke = -l.OdO 
flux = fluxr 
call shoot (pexit) 
b = choke 
choke = -l.OdO 
test = a*b 

if (test.lt. (O.dO) ) go to 24 


fluxr = fluxl 
fluxl = temp 

fluxl = (fluxl + fluxr) /2.0d0 
go to 23 

24 check = dab s (fluxl - fluxr) /fluxl 
cflux = fluxl 
pchoke = pexit 
if (check. It .bound) go to 25 
tenp = fluxl 

fluxl = (fluxl + fluxr) /2.0d0 
go to 23 
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25 


continue 


c Test now whether the back pressure desired can be matched. If 

c it can not be matched store the profiles for the choked flow, If 

c it can be matched proceed to find the proper mass flux and 

c exit pressure . 

if (pback. ge.pchoke) go to 26 

call output (iplot, lease) 
go to 100 

26 continue 

c Input the first two guesses for the mass flux. 

fluxOO = floor 
fluxO = cf lux 

c Calculate the exit pressures for the guessed mass flux. Abort 

c the procedure if the initial guess for flux gives choked flow. 

flux = fluxOO 

call shoot (pendOO) 

if (choke. ne. (1 .d0) ) go to 30 

write (*, 9) 

return 


30 


cc 


flux = fluxO 
iterat = iterat + 1 


cc 

call shoot (pendO) 


c Correct the guess for flux. If it was found that the last guess 

c for flux gave an infinite or negative pressure gradient in the 
c computational domain, average that guess with the previous one, 
c reset the choked flow flag and shoot again. 


if (choke. eq. (-l.dO) ) go to 40 
fluxO = 0 . 5d0* (f luxO+f luxOO) 
choke = -l.OdO 
go to 30 

40 deltaf = (pback - pendO) * (fluxO - fluxOO) / (pendO - pendOO) 
fluxl = fluxO + deltaf 
fluxOO = fluxO 
fluxO = fluxl 
pendOO = pendO 

c test for convergence. 

check = dabs (pendO -pback) /pback 
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if (iterat .ge.100) then 

write (*, ' ' program did not converge, returned to caller. ) ) 

return 

end if 

if (check. ge. error) go to 30 

c call output routine to store information for printing and 
plotting. 

cc write (*,'(" iterations = " , i5) * ) iterat 
call output (iplot, icase) 

100 hrin = hrin/hstep 
continue 

c 

c end of the big do loop. 


c Format statements . 

9 format (1 lx, 'Note: Initial guess for flux gives choked flow') 

60 return 
end 


c 


subroutine dsat<3? (dvfdp, dvfgc^p, difcfc>, difg<$>, press) 


c This subroutine calculates derivatives of certain saturation 

c properties as functions of pressure. The data is stored in the 
c array ddatdp(i, j,k) where i denotes the fluid property as follows: 

c l:d(vf)/dp 2:d(vfg)/dp 3:d(if)/dp 4:d(ifg)/dp 

c 

c For a particular j, k=2,3,4 contains the coefficients of a 

c quadratic which is the derivative of the cubic spline for the 

c curve of property type i (e.g. if) between point j and j+1. 
c 

c 


implicit real*8 (a-h,o-z) 

common /gas / gamma, cpgas, tcrit, vcrit, per it 

common /geomry/ hmean, radius, rin, rout, r inlet, dhdr, nspace, 

$ delta, epsi, pi, rmean 

common /flow / flux, omega, pinfty, pback, iinfty, tinfty, choke, 

$ cf lux, pchoke, closs, signum, ifluid, psmall 

common /thermo/ satdat (11,50, 4) , ndata, jprops, ddatdp(5, 60,3) , 

$ visdat (65,15), pvisc(15), tvisc(65), irows, jcols, 

$ tmin, tmax, pnin, pmax 
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common /var / p(120), h(120) 
common /stress/ tauzr, tauzth, v, mu 


real*8 mu, iinfty 

c Jump to the calculations if pressure is still between p(j) and 

P(j+1) . 

if ((press .ge. satdat (6, jprops, 1) ) .and. 

$ (press .le. satdat (6, jprops+1, 1) ) ) go to 200 

c Find the correct j for this property. 

jmove = int (dsign (1 . OdO, press - satdat (6, jprops, 1) ) ) 

100 jprops = jprops + jmove 

if ((press .It. satdat (6, jprops, 1) ) .or. 

$ (press .gt. satdat (6, jprcps+1, 1) ) ) go to 100 

c Calculate the property from the cubic spline derivative. 


200 z = (press - satdat (6, jprops, 1) ) 
dvfc£> = ddatcip (1, jprops, 1) + z * 

$ + z * 

dvfgc£> = ddatcfc (2, jprops, 1) + z * 

$ + z * 

difdp = ddatdp (3, jprops, 1) + z * 

$ + z * 

difgdp = ddatdp (4, jprops, 1) + z * 

$ + z * 


(ddatdp (1, jprops, 2) 
(ddatdp (1, jprops, 3) ) ) 
(ddatcip (2, jprops, 2) 
(ddatcip (2, jprops, 3) ) ) 
(ddatdp (3, jprops, 2) 
(ddatdp (3, jprops, 3) ) ) 
(ddatcip (4, jprops, 2 ) 
(ddatoip (4, jprops, 3) ) ) 


return 

end 


c 

c subroutine grad 

c 

c 

c This subroutine finds fluid prcperities and their derivatives 

c with pressure and calculates the axial gradients of pressure and 

c enthalpy derived from the equations of notion and energy, 

c Multiple entry points are used here to calculate the appropriate 

c properties for three separate flow regimes: sub-cooled liquid, 
c saturated liquid - vapor mixture, and superheated vapor, 
cc 

cc Multiple entry points are, unfortunately, not supported by 

cc Microsoft FORTRAN compiler version 3.31. The three entry points 

cc in the original subroutine are broken up into three separate 

cc subroutines, 

cc - Stephen Lau 

c 

c 
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cc 


implicit real*8 (a-h,o-z) 

^v, cairnon /gas / garrma, cpgas, tcrit, vcrit, pcrit 

cc common /geomry/ hmean, radius, rin, rout, rinlet, dhdr, nspace, 

cc $ delta, epsi, pi, rmean 

cc cormon /flow / flux, omega, pinfty, pback, linfty, tinfty, choke, 
cc $ cf lux, pchoke, closs, signum, lfluid, psmall 

cc cannon /thermo/ satdat (11,50,4), ndata, jprcps, ddatdp(5, 60,3) , 
cc $ visdat (65, 15) , pvisc(15), tvisc(65), irows, jcols, 

cc $ tmin, tmax, pmin, pnax 

cc cannon /var / p(120), h(120) 
cc cannon /stress/ tauzr, tauzth, v, mu 

cc real*8 mu, iinfty 

c 


subroutine liquid (press, enthpy,dpdr, didr) 


implicit real*8 (a-h, o-z) 

cannon /gas / gamna, cpgas, tcrit, vcrit, pcrit 

cannon /geomry/ hmean, radius, rin, rout, rinlet, dhdr, nspace, 

$ delta, epsi, pi, rmean 

cannon /flow / flux, omega, pinfty, pback, iinfty, tinfty, choke, 
$ cflux, pchoke, closs, signum, ifluid, psmall 

cannon /thermo/ satdat (11, 50, 4) , ndata, jprcps, ddatdp(5, 60,3) , 

$ visdat (65,15) , pvisc(15), tvisc(65), irows, ijcols, 

$ tmin, tmax, pcnin, pnax 

cannon /var / p(120), h(120) 
common /stress/ tauzr, tauzth, v, mu 

real*8 mu, iinfty 

Find the local film thickness. 

hfilm = hmean + (radius - rmean) *dhdr 

Find the mass velocity, g. 

g = signum* flux/ (2. OdO*pi*radius*hfilm) 

Take the viscosity and specific volume of the sub-cooled liquid 
to be that of the saturated liquid at the same temperature. 

temp = sat (ll,enthpy) 
v = sat (10, temp) 
mu = sat (7, temp) 
call shear 

thing = g*g*v/ radius* (l.OdO + radius/hfilm*dhdr) 
terml = - 2 . 0 d 0 *signum/hfilm*tauzr 
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term2 = cmega*cmega*radius/ (3.0d0*v) 
term3 = thing 


<3pdr = terml + term2 + term3 

terml = radius*cmega/ (g*hf ilm) *tauzth 
term2 = -omega * omega *radius/3.0d0 
term3 = v*thing 


didr = terml + term2 + term3 

return 

end 


c 


subroutine phas2t (press , enthpy , dpdr , didr) 


implicit real*8 (a-h, o-z) 

cannon /gas / gamma, cpgas, ter it, verit, perit 

common /geomry/ hmean, radius, rin, rout, runlet, dhdr, nspace, 

$ delta, epsi, pi, mean 

cannon /flow / flux, omega, pinfty, {back, iinfty, tinfty, choke, 
$ cflux, pchoke, closs, signum, lfluid, psmall 

cannon /thermo/ satdat (11,50,4) , ndata, jprpps, ddatdp(5, 60, 3) , 

$ visdat (65,15) , pvisc(15), tvisc(65), irows, ^cols, 

$ tmin, tmax, pmin, pnax 

cannon /var / p(120), h(120) 
cannon /stress/ tauzr, tauzth, v, mu 

real*8 mu, iinfty 

o Find the local film thickness. 

hfilm = hmean + (radius - mean) *dhdr 
c Find the mass velocity, g. 

g = signum* flux/ ( 2 . 0 d 0 *pi*radius*hfilm) 
c Find the saturation temperature, 

temp = sat (1, press) 

c Find vfg, if, and ifg for the given pressure. 

vfg = sat (3, press) 
if = sat (4, press) 
ifg = sat (5, press) 

c Calculate the quality. 
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qualty = (enthpy - if ) /ifg 

c Find the specific voluire of the saturated mixture at a given 

c pressure . 

vf = sat (2, press) 
v = vf + qualty*vfg 

c Find the derivatives of vf, vfg, if, and ifg with pressure. 

call dsatdp( dvfdp, dvfgdp, difdp, difgdp, press) 

c Find the derivatives of volume and enthalpy with respect to 

c pressure at constant quality . 

dvdp = dvfdp + qualty*dvfg<^> 
dicfc = difdp + qualty*difgc$) 

c Find the mixture viscosity by weighting the liquid and vapor 

c viscosities with the volume fraction. 

2 = ^(qualty*vg*sat (8,terrp) + (l.dO - qualty) *vf*sat (7, temp) ) /v 

c Calculate the denominator of the gradient expressions and 

c check for a singularity. 

phi = -g*g* (dvdp + vfg/ifg* (v - didp) ) 

if (phi.lt. (l.dO)) go to 100 

choke = l.OdO 

return 

100 continue 

c Find the shear stresses, 

call shear 

c Calculate the pressure and enthalpy gradients. 

quant = l.OdO + radius/hfilm*dhdr 
factor = g*g*v/radius 
quodl = g*g*vfg/ifg 

quod2 = -omega*radius*tauzth/ (g*hfilm) 
quod3 = omega* omega* radius / 1 . 5d0 
quod4 = - 2 . 0 d 0 *signum*v*tauzr/hfilm 

terml = factor*quant 
term2 = onega*OTiega*radius/ (3.0d0*v) 
term3 = - 2 . 0 d 0 *signum*tauzr/hfilm 
term4 = quodl* (quod2+quod3+quod4) 
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dpdr = ( t erml +t erm2 +term3+term4 ) / (1 . OdO - phi) 


factor = -l.OdO/ifg 

terml = -cmega*radius*tauzth/ (g*hfilm) 

term2 = cirega*omega*radius/l . 5d0 

term3 = - (v - Hi rip) *dpdr 

term4 = -2 . OdO*signum*v*tauzr/hf ilm 

dqdr = factor* (terml +term2+term3+term4) 

didr = (difclp + qua 1 ty *dif gdp) *dpdr + ifg*dqdr 

return 

end 


c 


subroutine vapor (press, enthpy, dpdr, didr) 


c 


implicit real*8 (a-h,o-z) 

common /gas / gamma, cpgas, ter it, verit, per it 

common /geomry/ hmean, radius, rin, rout, rinlet, dhdr, nspace, 

$ delta, epsi, pi, rmean 

common /flow / flux, omega, pinfty, pback, iinfty, tinfty, choke, 
$ cflux, pchoke, closs, signum, ifluid, psmall 

common /thermo/ satdat (11,50, 4) , ndata, jpreps, ddat<3p(5,60,3) , 

$ visdat(65, 15), pvisc(15), tvisc(65), irows, jcols, 

$ tmin, tmax, pmin, pmax 

common /var / p(120), h(120) 
common /stress/ tauzr, tauzth, v, mu 

real *8 mu, iinfty 

c Find the local film thickness. 

hfilm = hmean + (radius - rmean) *dhdr 
c Find the mass velocity, g. 

g = signum* flux/ (2.0d0*pi*radius*hfilm) 
c Find the temperature and determine the viscosity, 

temp = sat (1, press) 

$ + (enthpy - sat (4, press) - sat (5, press) ) /cpgas 

mu = superv (press, temp) 

c Find the specific volume (through ideal gas law) . 

r = cpgas* (gamma - 1. OdO) /gamma 
v = r*temp/press 

c Determine the mach number and check for a singularity. 
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machsq = g*g*v/ (gaitma*press) 
if (machsq.lt. (l.dO) ) go to 200 
choke = l.OdO 
return 


200 continue 

Find the shear stresses, 
call shear 


Calculate the pressure and enthalpy gradients. 


terml = 

$ 

term2 = 
term3 = 

$ 

term4 = 
<3pdr = 

terml = 
term2 = 
term3 = 
term4 
didr - 


-2 . OdO*signum*tauzr/hf ilm* (1 . OdO + 

(qarrma - l.OdO) *machsq) .. , 

: -(gaitma - l.OdO) * radius* anega*machsq*tauzth/ (g hfilm v) 

: anega*omega*radius/(3.0d0*v) 

* (l.OdO + 2 . OdO* (gamma - l.OdO) *machsq) 

= gamma*press*machsq/ radius* (1 . OdO + radius/hfilm dhdr) 

. 4 term2 + ?erm3 + term4)/(1.0d0 - machsq) 

= radius*omega/ (g*hf ilm) *tauzth 

= omega* omega* radius/ 1 .5d0 ^ 

= - 2 . 0 d 0 *signum*v*tauzr/hfilm 

= (gamma + l.OdO) / (gamma - l.OdO) *v*3pdr 
terml + term2 + term3 + term4 


return 

end 


c 


subroutine output (iplot, icase) 


c 

c 

c 

c 

c 

c 

c 

c 

c 

cc 

cc 

cc 

cc 

cc 

cc 

cc 

c 

c 


This subroutine prints out point by ^ information about 
temperature, pressure, and quality in the file face.dat . 

< fTS Sriaible iplot = 1 , this routine also generates 
a mi ^£tf!da£ Which be further processed by the program 
tof fom Ms post-processor prepares a plot of pressure, 
temperature, and quality profiles on a single graph. 

Most of the output statements in the routine are consented out. 

S^rS^i^^oS'^rr^nSiS qualS/at the 

"&&£££ ^^o^r^|« 0 sta^ts 

for generating the plot file are also commented out. 

- Stephen Lau 
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implicit real*8 (a-h, o-z) 

common /gas / gamma, cpgas, tcrit, vcrit, pcrit 

common /geomry/ hmean, radius, rin, rout, r inlet, dhdr, nspace, 

$ delta, epsi, pi, rmean 

common /flow / flux, omega, pinfty, pback, iinfty, tinfty, choke, 

$ cflux, pchoke, closs, signum, ifluid, psmall 

conmon /thermo/ satdat (11,50,4) , ndata, jprcps, ddatdp(5, 60, 3) , 

$ visdat (65,15), pvisc(15), tvisc(65), irows, jcols, 

$ tmin, tmax, pmin, pmax 

common /var / p(120), h(120) 
common /stress/ tauzr, tauzth, v, mu 

real*8 mu, iinfty, leak 

dimension r(120), t(120), qualty(l20) 

c Back out properties of interest from the enthalpy array. 

do 10 i=l,nspace+l 

check = (h(i) - sat (4,p(i) ) ) /sat (5,p(i) ) 
if (check.lt. (O.dO)) check = -1 
check = dint (check) 
if (check) 11, 12, 13 

11 qualty(i) = O.OdO 

t (i) = sat (ll,h(i) ) 
go to 10 

12 qualty(i) = (h(i) - sat (4,p(i) ) ) /sat (5,p(i) ) 
t (i) = sat (l,p(i) ) 

go to 10 

13 qualty(i) = l.OdO 

t(i) = sat (l,p(i) ) + (h(i)-sat (5,p(i) )-sat (4,p(i) ) ) /cpgas 
10 continue 

c Calculate the leakage rate and maximum leakage rate. 

leak = flux 
cleak = cflux 

c calculate shaft rpm. 

rpm = omega*30 . OdO/pi 

c Report if the flow was choked. 

iflow = 0 

if (pback. le. pchoke) iflow = 1 

c Integrate the pressure over the seal area by the trapezoidal 

c rule to find the opening force 

delr = delta 
r(l) = rinlet 
sum = O.OdO 
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do 31 i=2,nspace+l 
r (i) = r(i-l) + delr 
tleft = p(i-l) *r (i-1) 
tright = p(i)*r(i) 

sum = sum + signum*delr*pi* (tleft + tright) 

31 delr = delr*epsi 

force = sum 

c Check circumferential and radial flow reynold's numbers, 

c Print out the minimum values of each for reference. 

reyc = 1.0d+20 
reyr = 1.0d+20 

delr = delta 

do 1000 i = l,nspace+l 

radius = rinlet + delr*float (i-1) 

hfilm = hmean + (radius - rmean) *dhdr 

g = flux/ (2.0d0*pi*radius*hfilm) 

vf = sat (2,p(i) ) 

vfg = sat (3,p(i) ) 

v = vf + qualty (i) *vfg 

den = l.OdO/v 

vg = vfg + vf 

vise = (qualty (i) *vg*sat (8, t (i) ) 

$ + (l.OdO - qualty (i))*vf* sat (7, t(i)))/v 

rr = 2.0d0*g*hfilm/visc 

rc = den*omega*radius*hfilm/ (2 . 0d0*visc) 
if (rr.lt .reyr) reyr = rr 
if (rc.lt. reyc) reyc = rc 
delr = delr*epsi 
1000 continue 

c Store output information in files for printing and plotting. 

cc open (unit =15, file=* test .dat ' , access=' sequential ' , status='new' ) 

cc write (15, 20) 

cc write (15, 21) 

cc if (ifluid.eq.l) write (15, 36) 

cc if (ifluid.eq.2) write (15, 46) 

cc if (ifluid.eq.3) write (15, 47) 

cc if (ifluid.eq.4) write (15, 48) 

cc write (15, 37) icase 

cc write (15, 22) rin 

cc write (15, 23) rout 

cc write (15, 24) hmean 

cc write (15, 42) dhdr 

cc write (15, 25) rptn 

cc write (15, 26) pinfty 

cc write (15, 27) pback 
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cc write (15, 28) tinfty 

cc write (15, 38) closs 

cc write (15, 29) leak 

cc write (15, 33) cleak 

cc write (15, 34) pchoke 

cc write (15, 45) reyr 

cc write (15, 44) reyc 

cc write (15, 43) force 

cc if (iflow.eq. (l.dO) ) write(15,39) 

cc write (15, 30) 

cc do 35 i=l,nspace+l 

cc 35 write (15, 32) i, r(i), t(i), p(i), qualtyd) 

cc write (15, *(///)') 

cc close (unit=15) 

cc if (iplot.eq.O) go to 40 

cc open (unit=14 , f ile= ' fprof.dat ' , access= ' sequential ’ , status= ' old ' ) 

cc write (14,*) icase 

cc write (14,*) ifluid 

cc write (14,*) rin 

cc write (14,*) rout 

cc write (14,*) hmean 

cc write (14,*) dhdr 

cc write (14,*) rpm 

cc write (14,*) pinfty 

cc write (14,*) pback 

cc write (14,*) tinfty 

cc write (14,*) leak 

cc write (14,*) nspace 

cc write (14,*) closs 

cc write (14,*) force 

cc write (14,*) iflow 

cc do 41 i = l,nspace+l 

cc write (14,*) r(i) 

cc write (14,*) t(i) 

cc write (14,*) p(i) 

cc 41 write (14,*) qualty(i) 

cc close (unit=14) 

c+ " 

c Tabulate results in format similar to Iau's laminar program, 
c 

rmean = 0.5d0*(rin + rout) 

hrin = hmean - 0 . 5d0*dhdr* (rout - rin) 

write (12,501) hrin, reyc, reyr, qualty (nspace+1) , leak*signum, force 
501 format (3 (lx, lpel2 .3) , lx, OpflO. 3, 2 (lx, lpe!2.3) ) 



c 

c+ 


c 

■+c 


40 return 


c Format statements- 

cc 20 format (1 Ox, lh 'Steady State Analysis for a Face Shaft 

« y&grr**'''*** ,in 


cc 21 format (5x, 


/) 


cc 22 format (5x, 
cc 23 format (5x, 1 
cc 24 format (5x, ' 
cc 25 format (5x, ' 
cc 26 format (5x, ' 
cc 27 format (5x, ' 
cc 28 format (5x, ' 
cc $ ’ 

cc 29 format (5x, ’ 
cc $ 

cc 30 format (lx, 
cc $ 


Inner Seal Radius = ' , lpel0.3, ' meters ) 
Outer Seal Radius = 1 , IpelO . 3, ' meters ) 
Mean Film Thickness = ',lpel0.3, - meters') 
Shaft Rotation Speed = IpelO. 3,' fF™' ) _ 
Reservoir Pressure = ', IpelO. 3, Pascals ) 
Back Pressure = IpelO. 3,' Pascals ) 
Reservoir Temperature = *, IpelO. 3, 

degrees Kelvin') 

Leakage Rate — , lpel2 .5, 

kg/ second') 

Point R . Temperature 

Pressure Quality ' / ) 


cc 32 format (2x,i3,7x,lfel0. 4 7x lpel0.4 


' 1 ' , . io c. 

Maximum Leakage Rate - , lpelz . o, 

kg/second') . _ 

Exit Pressure for Choked Flow - , IpelO. 3, 

Pascals ' ) Sealed Fluid is WATER' /) 

case # \i3/) 

Inlet Loss Coefficient = ’,el0.3) 

Flow is Choked'/) 

Coning Angle = ', IpelO. 3,' radians') 
Opening Force = \lpel2.5,’ Newtons’/) 
~ SSSSS 'Minimum CiJSial Re - ;.lpel0.3> 
cc 45 format (5x, ' .W'phRh - 

CC S x ' Sealed Fluid is LOX - GCK-/> 

cc 47 format lx, _ sealed Fluid is NITROGEN’/) 

cc 48 format (lx, 

end 


cc 33 format (5x, 
cc $ 

cc 34 format (lx, 
cc $ 

cc 36 format (lx, 
cc 37 format (lx, 
cc 38 format (5x, 
cc 39 format (5x, 
cc 42 format (5x, 
cc 43 format (5x, 


c 


subroutine runge (n,y,f,r,h,m,k) 


c — 

c 

c 

c 

c 

c 

c 

c — 


this subroutine performs numerical quadrature of coupled dif 
f a r^nt i a 1 eauations by gill's method, a variant of the Runge-Kutta 
technique.^This routine was drawn from appendix c of "Viscous 
Fluid Flow” by Frank M. White. 


implicit real *8 (a-h,o z) 
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dimension y(10), f(10), q(10) 
m = m + 1 

go to (1, 4,5, 3,7), m 

1 do 2 i = l,n 

2 q(i) = O.OdO 
a = 0 . 5d0 
go to 9 

3 a = l.OdO + dsqrt (0.5d0) 

4 r = r + 0.5d0*h 

5 do 6 i = l,n 

v(i) = v(i) + a*(f(i)*h - qd)) „ _ ,. v 

6 q(i) = 2.0d0*a*h*f (i) + (l.OdO - 3.0d0*a) q(i) 
a = l.OdO - dsqrt (0.5d0) 

go to 9 

7 do 8 i — ^ 

8 y (i) = y (i) + h*f(i)/6.0d0 - q(i)/3.0d0 
m = 0 

k = 2 
go to 10 

9 k = 1 
10 return 

end 


c 


subroutine shear 



This subroutine corputes the radial (TftUZR) and circumferential 
(TAUZTH) shear stresses at the seal faces- Expressions or 
both shear stresses are in a form: 


- - Cl In (P) - C2 = 0 

P 

where p is the square root of the friction factor. The 
stresses may be recovered from the formula: 


Twall = (PU) / (8v) 

here U is the mean velocity and v is the specific volume . 

-he constants cl and c2 are determined by the fluid properties, 
he ^raSng conditions of the seal, and the stress component 
aider consideration. Solution of this transcendental 
equation is accomplished by Newton’s method- 


implicit real* 8 (a-h,o-z) 
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common 

common 

$ 

cannon 

$ 

cannon 

$ 

$ 

cannon 

common 


/gas / 
/geomry/ 

/flow / 

/thermo/ 


/var / 
/stress/ 


KS nspace, 

iirtty, t ****<*•. 

cf lux, pchoke, closs, signum, if ^ x f' ... 

satdat (11,50,4) , ndata, Dprops, ddatdp (5^6°, 3) , 
visdat (65, 15) , pvisc(15), tvisc(65), irows, oco , 
tmin, tmax, pmin, pnax 
p (120) , h(120) 
tauzr. tauzth, v, mu 


real*8 mu, iinfty 

c Find the local film thickness. 

hfilm = hmean + (radius - rmean) *dhdr 


q Find the mass velocity, g. 

g = flux/ ( 2 . 0 d 0 *pi*radius*hfilm) 

c set error bound for iteration. 

error = 1.0d - 06 


c Solve for axial shear stress. 


tauzr = dabs (tauzr) 
arg = 8.0d0*tauzr/ (g*g*v) 
prs = dsqrt(arg) 
cl = 0. 8838d0 

c 2 = 0.142d0 + cl*dlog (g*hf ilm/rnu) 
pcrit = dexp (-c2/cl) . . in 

if (prs .It. pcrit) prs = pent lO.OdO 

100 _ 2.0d0 + cl*pold* (l.OdO - dlog(pold)) - c2*pold 

bot = l.OdO/pold + cl 

i?(prs t a^ b °^crit) prs = dsqrt (pold*pcrit) 
if (dabs (pold - prs) / prs .gt. error) go to 

Limit the friction factor for transition flows. 


reyn = g*hfilm/mu 
if (reyn.lt. (4. 0d3) ) prs = 

temp = prs*g 

tauzr = temp*temp*v/8.0d0 


0.1694d0 


c Solve for circumferential shear stress. 


if (omega .eg. (O.dO)) return 
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tauzth = dabs (tauzth) 

u = radius* omega 

arg = 8. OdO*v* tauzth/ (u*u) 

prs = dsqrt (arg) 

cl = 1.763d0 

arg = u*hfilm/ (mu*v) 

c2 = 0.83d0 + cl*dlog(arg) 

per it = dexp(-c2/cl) 

if (prs .It. perit) prs = pent * lO.OdO 


.OdO - dlog(pold)) - c2*pold 


200 pold = prs 

top = 2. OdO + cl*pold* (1 . 
bot = 1. OdO /pold + cl 
prs = top/bot 

if (prs .It. perit) prs = dsqrt (pold*pcnt) 
if (dabs (pold - prs) / prs .gt. error) go to 200 


temp = prs*u 

tauzth = temp*tenp/ (8.0d0*v) 


return 

end 


e 


subroutine shoot (pexit) 


c — 

c 

c 

c 

c 

c 

c — 


c 

c 


This function returns the value of the pressure at the exit 
boundary of the computational domain and records the pressure 
and enthalpy at each finite difference point in that domain. 


implicit real*8 (a-h, o-z) 

common /gas / gamma, qpgas, tcrit, verit, perit 
common /geomry/ hmean, radius, run, rout, r inlet, dhdr, nspace, 
c delta, epsi, pi, rmean 

cannon /flow / flux, anega, pinfty, ^^^il ' 

$ cf lux, pchoke, closs, signum, psmall 

cannon /thermo/ satdat(U,50,4>, nda ’;^, icils 

§ visdat (65, 15) , pvisc(15), tvisc(65), irows, ^cois, 

$ tmin, tmax, pmin, pmax 

common /var / p(120), h(120) 
conmon /stress/ tauzr, tauzth, v, mu 

real*8 mu, iinfty, iinlet 

dimension y (10) , f(10) 

Define number of equations to be solved and initialize para- 
meters for integrat ion scheme . 
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ra = 0 
n = 2 

radius = rinlet 

c Find the local film thickness. 

hfilm = hmean + (radius - rmean) *dhdr 
c Find the mass velocity, g. 

g = signum*flux/( 2 . 0 d 0 *pi*radius*hfilm) 

c Calculate the entrance losses and modify inlet pressure and 

c enthalpy. 

c Define error bound for iteration. 

error = 1.0d-04 

c Make initial guess that the end state corresponds to pure liquid. 

volume = sat (10,tinfty) 
thing = g*volume 

-i-inlpt = iinftv - thing* thing/2. OdO 

pinlet = pinfty - (closs + l.OdO) /volume*thing*thing/2 . OdO 

c Check for unrealistic pressure loss at the inlet. 

if (pinlet. ge.psmall) go to 3 
choke = l.OdO 
return 
3 continue 

c Check if end state is indeed pure liquid. 

test = sat (4, pinlet) 
if (iinlet.le.test) go to 5 

c Iterate for saturated mixture end state. 

voll = volume 


1 ptemp = pinlet 

volf = sat (2, pinlet) 
volfg = sat (3, pinlet) 
hf = sat (4, pinlet) 

hfg = sat (5, pinlet) 
qualty = (i inlet - hf ) /hfg 
vol2 = volf + qualty*volfg 


thing = g*vol2 

i in let = iinfty - thing*thing/2 . OdO 
vol = (voll + vol2)/2.0d0 
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pinlet = pinfty - (closs + 1 . OdO) *thing*thing/ (2 . OdO*vol) 

Test for unrealistic pressure loss at inlet. 

if (pinlet. ge.psmall) go to 2 

choke = l.OdO 

return 

Check for convergence. 

2 test = dabs (pinlet -ptenp) /pinlet 
if (test. ge. error) go to 1 

5 continue 

Initialize index to mark variable arrays and apply the initial 
conditions in the shooting algorithm. 

jpoint = 1 
delr = delta 
y (1) = pinlet 
y (2) = iinlet 

Store the current values of pressure and enthalpy and integrate 
both the axial momentum and energy equations. 

8 p (jpoint) = y(l) 
h (jpoint) = y(2) 

if (jpoint. eq. (nspace+1) ) go to 7 
6 call runge (n,y, f, radius, delr, m,k) 
press = y (1) 

Abort procedure if pressure goes below known saturation data. 

if (press. ge.psmall) go to 4 
choke = l.OdO 
go to 18 

4 enthpy = y(2) 
goto (10,20), k 

10 check = (y(2) - sat (4,y (1) ) ) /sat (5,y (1) ) 
if (check.lt. (O.dO) ) check = -l.OdO 
check = dint (check) 
if (check) 12, 14, 16 
12 call liquid (press, enthpy, dpdr, didr) 
f(l) = dpdr 
f (2) = didr 

go to 6 , . 

14 call phas2t (press, enthpy, c^pdr, didr) 
if (choke. eq. (l.dO) ) go to 18 
f(l) = dpdr 
f (2) = didr 

go to 6 . . 

16 call vapor (press, enthpy, dpdr, didr) 



if (choke. eq. (l.dO) ) go to 18 
f (1) = dpdr 
f (2) = didr 
go to 6 

20 jpoint = jpoint + 1 
delr = epsi*delr 
go to 8 

Give pexit the pressure calculated at the end of the ccrrpu 
tational domain. 

7 pexit = y (1) 

18 return 
end 


real *8 function superv (press, temp) 


c 

c this function computes superheated vapor viscosity by interpol- 
c ating from from a table. 

c 

c 


implicit real*8 (a-h, o z) 

common /gas / 
common /geomry/ 

$ 

common /flow / 

$ 

common /thermo/ 


$ 

$ 


common /var / 
nrrrmnn /stress/ 


SE SSus. tC rS'rSt1 t riSS' <Mr. nspace, 

SS?'oSa: P^fSTSaack, iinfty, tinfty, choke, 
cf lux, pchoke, closs, signum, iflurd, psmali 
satdat (11,50,4) , ndata, Dprops, ddatcp(5, 60, J) , 
visdat (65^ 15) , pvisc (15) , tvisc(65), irows, ^cols, 
tmin, tmax, pmin, pmax 
p(120) , h(120) 
tauzr, tauzth, v, mu 


c 

c 


c 


real* 8 mu, iinfty 

Test to determine whether thermodynamic state is within the 
bounds of the table. 


serror = O.OdO 
if ( (press . It . pmin) . or 
if ((terrp.lt. tmin) .or. (terp.ge 
if (serror .eq. (1 .d0) ) go to 100 


(press .ge. pmax)) serror = l.OdO 
tmax)) serror = l.OdO 


Find the appropriate pressure range. 


k = 0 

10 k = k + 1 
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if (press. ge. pvisc (k) ) go to 10 
j = k - 1 

c Find the appropriate temperature range. 

k = 0 

20 k = k + 1 

if (temp.ge.tvisc (k) ) go to 20 
i = k - 1 

c — Test for proximity to the saturation dome. 

if (tenp.ge.tcrit) go to 40 

psat = sat (6, temp) 

if (pvisc ( j+1) .gt.psat) go to 30 

c — Bilinear interpolation by inverse-area rule. 

40 al = terrp - tvisc(i) 
a2 = tvisc (i+1) - tenp 
bl = press - pvisc (j) 
b2 = pvisc (j+1) - press 

factl = al*bl 
fact2 = a2*bl 
fact 3 = a2*b2 
fact 4 = al*b2 

factor = 1.0d0/ (factl + fact2 + fact3 + fact4) 

superv = factor* (fact3*visdat (if j) + fact4*visdat (i+lf j) 
$ + factl*visdat(i+l, j+1) + fact2*visdat (i, j+1) ) 


go to 200 

c Special interpolation routine near the saturation done. 

30 vhigh = sat (8, temp) 
tempi = sat (1, pvisc (j) ) 
vlowl = sat (8, tempi) 
vlow2 = visdat (i+1/ j) 
factl = temp - tempi 
fact2 = tvisc (i+1) - temp 
factor = 1.0d0/ (factl + fact2) 
vlow = factor* (fact2*vlowl + factl*vlow2) 

factl = press - pvisc (j) 
fact2 = psat - press 
factor = 1.0d0/ (factl + fact2) 

1000 continue 

superv = factor* (fact2*vlow + factl*vhigh) 
go to 200 
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100 write (5, 101) 

200 return 

c Format statement. 

101 format (lh 'Pt. of interest is out of bounds’) 

end 
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